SUBMITTED

1

Space Time MUSIC: Consistent Signal Subspace
Estimation for Wide-band Sensor Arrays

arXiv:1711.01631v2 [stat.AP] 1 Mar 2018

Elio D. Di Claudio, Raffaele Parisi, and Giovanni Jacovitti

Abstract—Wide-band Direction of Arrival (DOA) estimation
with sensor arrays is an essential task in sonar, radar, acoustics,
biomedical and multimedia applications. Many state of the
art wide-band DOA estimators coherently process frequency
binned array outputs by approximate Maximum Likelihood,
Weighted Subspace Fitting or focusing techniques. This paper
shows that bin signals obtained by filter-bank approaches do
not obey the finite rank narrow-band array model, because
spectral leakage and the change of the array response with
frequency within the bin create ghost sources dependent on the
particular realization of the source process. Therefore, existing
DOA estimators based on binning cannot claim consistency even
with the perfect knowledge of the array response. In this work,
a more realistic array model with a finite length of the sensor
impulse responses is assumed, which still has finite rank under
a space-time formulation. It is shown that signal subspaces at
arbitrary frequencies can be consistently recovered under mild
conditions by applying MUSIC-type (ST-MUSIC) estimators to
the dominant eigenvectors of the wide-band space-time sensor
cross-correlation matrix. A novel Maximum Likelihood based
ST-MUSIC subspace estimate is developed in order to recover
consistency. The number of sources active at each frequency
are estimated by Information Theoretic Criteria. The sample
ST-MUSIC subspaces can be fed to any subspace fitting DOA
estimator at single or multiple frequencies. Simulations confirm
that the new technique clearly outperforms binning approaches
at sufficiently high signal to noise ratio, when model mismatches
exceed the noise floor.
Index Terms—Direction finding, wide-band sensor arrays,
space-time processing, signal subspace, coherent focusing, MUSIC, Weighted Subspace Fitting, WAVES, AIC, BIC.

I. I NTRODUCTION

W

IDE-BAND signal processing with sensor arrays is
a relevant field of research in widespread remote
sensing applications. In particular, Direction of Arrival (DOA)
estimation in sonar, multi-band radar, acoustic surveillance,
seismics, video conferencing and Ultra-Wide Band (UWB)
radio mostly involves signal bandwidths that largely violate
the narrow-band assumption.
Wide-band DOA estimation is usually tackled by four main
approaches:
• Frequency binning [1], performed by decomposing the
wide-band signals into a set of narrow-band signals
(bins), each tuned to a different frequency. The spatial
information carried by a set of narrow-band spatial
covariance matrices (SCMs) is combined into a unique
DOA estimate by either a multi-frequency extension of
Manuscript received March 2, 2018; revised Month, Year.
The authors are with the Department of Information Engineering,
Electronics and Telecommunications, University of Rome “La Sapienza”,
via Eudossiana 18, 00184 Rome, Italy. (e-mail: {elio.diclaudio, raffaele.parisi}@uniroma1.it, gjacov@infocom.uniroma1.it)

Maximum Likelihood (ML) [1]–[3], MUSIC [4] and
Weighted Subspace Fitting (WSF) [5], [6] narrow-band
DOA estimators, or a coherent focusing stage [7]–[13].
• Time Difference of Arrival (TDOA) estimation, where
wavefronts are reconstructed by fitting estimates of the
differential time delays among pairs of sensors [14].
• Wide-band adaptive steered beamforming, optimized by
a Minimum Variance (MV) [15] or ML criterion in the
frequency domain [16].
• Compressed Sensing (CS), where the sparsity of sources
is enforced by a penalty functional in a quantized DOA
domain [5], [17].
The resolving capability of TDOA approaches is limited
by sampling and windowing and by the presence of multiple
sources and colored noise [14].
Wide-band steered beamforming and CS furnish very robust
but inconsistent DOA estimates, respectively asymptotically
biased and resolution limited.
Since parametric narrow-band DOA estimation is consistent
and can cope with multi-source, correlated noise and arbitrary
array geometries, there is a strong interest in extending its use
to the wide-band case through frequency binning. However
the observation time must be sufficiently small to comply with
the assumption of stationarity of the sources and to provide an
adequate number of independent snapshots for the stable SCM
estimation within each bin. These contrasting requirements
enforce an unavoidable spread of the binning filters in the
frequency domain [18].
The problem is sketched in Fig. 1 for a uniform linear array
(ULA). The spectral support of each noiseless plane wave
source propagating at speed u and impinging from azimuth
θ, referred to broadside, is a Dirac wall based on the line
kx = −ω sin (θ)/u in the plane formed by the wavenumber
kx and the angular frequency ω, whose mass is modulated by
the source spectrum [19].
In binning approaches, best alignments of DOA estimates
from different bins are inferred [5], [8]. However, the DOA
ambiguity caused by variations of the source spectrum and the
steering vector [20] within each bin and the spectral leakage
among adjacent bins due to window side-lobes cannot be
recovered, since the narrowband array model does not account
for them. In particular:
• Spatially spread, uncalibrated ghost sources emerge from
noise background for increasing signal to noise ratio
(SNR) and/or sample size [12], [21], [22].
• Therefore, the finite rank structure of the narrow-band
array model [20] does not hold within each bin.
• The model becomes unidentifiable and the detection of
the source number [23] is an ill-posed problem.

2

SUBMITTED

Fig. 1. The spectral support of four far field noiseless and uncorrelated
sources impinging on an ULA in the half plane whose coordinates are the
the wavenumber kx and the angular frequency ω > 0 [19]. The wavenumber
(and therefore the DOA) uncertainties due to the finite bin width and to a
typical focusing process [8] are illustrated.

The model mismatch generates a signal-dependent DOA
bias.
• For random sources, an excess DOA variance appears at
high SNR [12] and may even hamper the convergence of
optimal DOA estimators [24], such as ML [3] and WSF
[5], [6], [25].
This mismatch undermines the consistency of the existing
narrow-band DOA estimators even in the absence of noise used
on frequency bins and is especially relevant for wide spatial
apertures or when wide-band, short-time acquisitions are required. It adds up to calibration and focusing mismatches, all
relevant at high SNR [24], [26].
Consistency guarantees that the limit DOA solution for
vanishing noise or increasing the sample size is unique and
satisfies the underlying physical model, regardless of signal
realizations. This is important for scientific experiments, array
calibration, target tracking and nulling static interference in
telecommunications.
Our goal is to obtain a consistent estimate of the signal
subspaces, asymptotically free from artifacts to replace the
SCM counterparts. With reference to Fig. 1, it is clear that
this task requires a joint parametric modeling in space and
time.
In [27] high resolution, narrow-band SCMs were estimated by suppressing the spectral leakage by a bank of MV,
distortion-less, delay and sum (DS) beamformers, tuned to
the same frequency. This approach involves cross-temporal
moments through the space-time array covariance matrix
(STCM), which only requires time sampling of sensor outputs
at a sufficient rate. The STCM exhibits fast statistical convergence and low bias even for short data records. However, the
resulting Capon SCM estimator of [27] suffers from the small
sample issues of MV beamforming [28] and its structure was
not analyzed for DOA estimation purposes.
In this work, starting from the finite rank property of the
wide-band source signature on the STCM, under the assumption of finite lengths of the source to sensor impulse responses
[27], an extension of the narrow-band DOA identifiability
conditions [3], [20] shows that consistent wide-band DOA
estimation is theoretically feasible by space time ML or WSF
•

approaches. However, such estimators require the knowledge
of the number and the impulse responses of all sources and
a very accurate coarse DOA estimation, which is at least unpractical.
In [22], a MUSIC [20] type condition allowed to asymptotically recover finite rank narrow-band signal subspaces
at arbitrary frequencies from the the wide-band dominant
eigenvectors of the STCM (or any other statistic sharing a
similar structure [29], [30]), rather than from filtered array
data. Earlier versions of this generalized Space-Time MUSIC
(ST-MUSIC) [22], [31] subspace estimator did not allow for
a sound statistical estimation of the number of sources and
had a slowly decaying asymptotic bias, which hampered the
numerical stability of the DOA estimates, especially at extreme
SNR and with small samples, without the use of heuristic
regularization [32].
Therefore, a novel ST-MUSIC subspace estimator, sharing
similarities with the Signal Subspace-MUSIC (SS-MUSIC)
[33], was derived from the ML formulation of an inversion
problem. The ML setting optimizes the statistical efficiency
and increases the robustness of the sample signal subspace
at extreme SNR, thus allowing the development of subspace
rank estimators, based on Akaike [34] and Bayesian [23], [35]
Information Theoretic Criteria (AIC and BIC respectively).
A perturbative analysis of the ML ST-MUSIC signal subspace estimate established its consistency and the optimal
weighting in the WSF sense [6]. Therefore, the sample ML
ST-MUSIC subspaces at different frequencies can replace the
inconsistent SCM counterparts in any narrow-band or wideband subspace based DOA estimator [5], [11]–[13] regardless
of the spatial coherency of sources.
Numerical simulations of the ML ST-MUSIC subspaces
with WSF based DOA estimation schemes [6], [12] show a
definite statistical performance improvement at high SNR of
the ML ST-MUSIC subspaces over the classical counterparts,
approaching the Cramer Rao DOA variance bound, at the cost
of a slightly increased low SNR threshold in some cases,
whose causes are discussed, and higher computation burden.
After a notation Sect. II, the finite length wide-band array
response model is analyzed in Sect. III for DOA estimation
purposes. In Sect. IV the ST-MUSIC concepts are developed.
In Sect. V, the ML ST-MUSIC subspace estimator is developed and AIC and BIC rank estimators are deduced, based on
the first order finite sample perturbation models of the STCM.
In Sect. VI the computational analysis is sketched. In Sect.
VII the validity of the ML ST-MUSIC model is assessed by
numerical simulations. Finally, conclusion is drawn in Sect.
VIII.
II. N OTATION
Throughout the paper matrices are indicated by boldface,
capital letters, vectors by boldface, lowercase letters. The
transpose of matrix A is indicated by AT , the conjugate by
A∗ , the Hermitian transpose by AH . The pseudo-inverse of A
is A† . IM is the square identity matrix of size M . The operator
diag{A} creates a column vector with the main diagonal of
A, diag{a} creates a diagonal matrix with the elements of

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

vector a placed on its main diagonal. trace(A) is the trace
of matrix A. det(A) is the determinant of A. A ⊗ B is
the Kronecker product between matrices A and B. δkl is the
Kronecker
√ symbol, equal to one if k = l and zero elsewhere.
j = −1 is the imaginary unit. Sub-matrices are indexed
by MATLAB-like conventions [36]. For instance, A(:, k) is
the k−th column of matrix A and A(1 : m, 1 : n) is the
upper left submatrix of A of size m × n. The operator ∗
indicates the temporal convolution of two signals. E {x} is the
expected value of the random variable x and var (x) denotes
its variance. The covariance between the random variables x
and y is indicated by cov (x, y) . The Frobenius norm [36] of
A is indicated by |A|F . Sample quantities are denoted by a
hat superscript (e.g., Â).
III. A RRAY MODEL
A sensor array with M sensors is immersed in a wavefield and receives the signals sd (t) (d = 1, 2, . . . , D), radiated
by D < M wide-band point sources, whose DOAs are
characterized by the generic coordinate vectors θd . Provided
that the sensors and the propagating medium are linear, the
signals xm (t) (m = 1, 2, . . . , M ), received by the M sensors
at time t, are modeled as
xm (t) =

D
X

hm (t; θd ) ∗ sd (t) + vm (t)

(1)

d=1

where hm (t; θd ) is the impulse response from the d-th source
to the m-th sensor, including the propagation channel (i.e.,
the Green’s function) and the front end filtering1 , and vm (t)
is the additive sensor noise term, assumed as statistically independent of source signals. All signal components can be either
complex envelopes with respect to an angular frequency ω0 , or
real valued base-band signals. For a down conversion carrier
angular frequency ω0 and a sampling period T , the array
response at the discrete time angular frequency ν coincides
with the response to a continuous time angular frequency
ω (ν) = ν/T + ω0 .
The sequences of consecutive P array output samples
xm (n) = xm (t)|t=nT for n = 1, 2, . . . , N are collected in a
space time snapshot (STS) of dimension M P × 1

T
xST (n) = xT1 (n) xT2 (n) · · · xTM (n)
(2)

3

common effective overall length Ld = max1≤m≤M Lm,d for
all hm (θd ) related to the generic d-th source.
Each impulse response hm (θd ), zero padded to Ld , is used
to build a Toeplitz convolution matrix of size P ×(P +Ld −1)


hm (θd )
0
···
0
 0
hm (θd ) · · ·
0 

(5)
Hm (θd ) = 
 ···
···
···
··· 
0
0
· · · hm (θd )
for m = 1, 2, . . . , M . Finally, the STS (2) is compactly
expressed as
xST (n) =

D
X

H (θd ) sd (n) + v (n)

(6)

d=1

where H (θd ) stacks Hm (θd ) for m = 1, · · · , M and is the
multi-channel convolution matrix of the sampled d-th source
signal

T
sd (n) = sd (n) sd (n − 1) · · · sd (n − P − Ld + 2)
(7)
and the M P × 1 vector

T
T
(n)
v (n) = v1T (n) v2T (n) · · · vM
(8)
collects the additive sensor noise samples

T
vm (n) = vm (n) vm (n − 1) · · · vm (n − P + 1)
.
(9)
Equations (5) and (6) show that the rank of H (θd ) is upper
bounded by P + Ld − 1 [27]. The SCM model [20] can be
recovered by assuming Ld = 1 and setting P = 1.
A. Space time covariance matrix
To tackle the consistency problem, cross-temporal moments are required to model the spectral spread. To this
purpose,
the Space-Time Covariance Matrix (STCM) RST =

[27] is a natural starting point for
E xST (n) xH
ST (n)
exploiting the model (6), because it collects all available
second-order information for zero mean, circularly complex
signals and noise with bounded fourth-order moments, if
P > max (Ld + Ls,d ), where Ls,d is the correlation length
d

of the dth source2 . It is worth to point out that in the same
scenario the SCM estimate would require the use of a DFT of
length Pw  P > max (Ld + Ls,d ), because its bias slowly
d

decays as Pw−1 [1], [8], [18].

T
In the sequel, the vector of source signals
xm (n) = xm (n) xm (n − 1) · · · xm (n − P + 1)
.

T
(3)
(10)
s (n) = sT1 (n) · · · sTD (n)
Each discrete time sensor impulse response hm (θd ) (m =
1, 2, . . . , M , d = 1, 2, . . . , D) represents the discrete time is assumed as a realization of a circular, wide sense stationary
process with zero mean and covariance
counterpart of hm (t; θd ) and is described by the row vector




S11 · · · S1D

hm (θd ) = hm,d (0) hm,d (1) · · · hm,d (Lm,d − 1) .
S = E s (n) sH (n) =  · · · · · · · · · 
(11)
(4)
SH
·
·
·
S
DD
1D
Its duration Lm,d may be infinite, if hm (θd ) contains poles

H
[38]. However, |hm,d (n)| typically quickly decays with n whose blocks Skl = E sk (n) sl (n) are the crosswell under the noise floor and it makes sense to assume a covariance matrices between each pair of impinging signals.
where

1 Following this formulation, h (t; θ ) may collect the effects of multiple
m
d
coherent reflections originated by the same driving signal sd (t) [37].

2 The use of other statistics based on (6), such as biased [22], structured
[39] or cyclic correlation matrices [29], is deferred to future work.

4

SUBMITTED

Sensor noise is assumed zero mean, circular
and inde
pendent of source signals, with STCM E v (n) vH (n) =
λv Rvv where λv ≥ 0. Then, the array STCM is

H
RST = E xST (n) xH
ST (n) = HST SHST + λv Rvv (12)


where HST = H (θ1 ) · · · H (θD ) .
B. Space time signal and noise subspaces
For known Rvv = Rv RH
v , the orthonormal and complementary bases Es for the wide-band signal subspace and Ev
for the wide-band noise subspace are classically defined by
the eigen-decomposition (EVD) [36] of the whitened STCM
−H
H
R−1
= Es Λs EH
v RST Rv
s + λv Ev Ev

(13)

where the diagonal matrix Λs = diag{λ1 , . . . , λη } contains
the dominant η eigenvalues in non-increasing order (i.e., λ1 ≥
λ2 ≥ . . . ≥ λη > λv ). In particular, the subspace Es has
dimension


D
X
η≤
(P + Ld − 1) ≤ D P + max (Ld ) − 1
(14)
d=1

d

and the column span of Es lies within the column span of
HST , so that [6]


C1



H (θ1 ) · · · H (θD )  ... 
Es = R−1
(15)
v
CD
where the blocks Cd for d = 1, 2, . . . , D have size
(P + Ld − 1)×η and point out the contributes of each source.
C. DOA identifiability from the STCM by ML and WSF
approaches
Model (12) is theoretically amenable to the straight extension of ML [1]–[3] and WSF [5], [6] narrow-band DOA
estimators3 , assuming for instance the exact knowledge of
the matrix impulse response H(θ) for all θ, of the noise
covariance Rvv up to a positive scalar, and of the exact number
of paths to search for. In addition, the source combination
which achieves the target RST must be unique [20].
In particular, the rank η of Es must remain smaller than
M P for the existence of Ev , crucial for the unambiguous identification of λv and DOAs [6], [20]. By (6), each uncorrelated
source contributes to η by a variable quantity ranging from one
(e.g., for a pure complex sinusoid) [22] up to P + Ld − 1 for
a random source spanning the full array bandwidth. Since for
D full-band sources η > DP , the number of resolvable wideband sources is strictly smaller than M . Further limitations on
the identifiability may be imposed by array ambiguities and
spatial coherence [20].
On the contrary, narrow-band sources generate only a few
significant eigenvalues and proper techniques [40]–[42] can
already estimate more DOAs than sensors if the source spectra
do not overlap.
3 For zero mean circular Gaussian signals and noise the STCM is a sufficient
statistic for DOA and S parameters of model (6).

However, the wide-band ML and WSF approaches based
on (6) lead to awkward, huge-dimensional problems that need
to be initialized by accurate (within half beam-width or less
[2]) prior DOA estimates for successful local convergence. In
particular, the difficulty of measuring the full H(θ) [38] motivated us to develop approaches based only on the knowledge
of the array harmonic response at a discrete set of temporal
frequencies.
IV. S PACE T IME MUSIC
In [22], a set of consistent and optimally weighted in
the WSF sense [6] narrow-band signal subspace bases at
arbitrary frequencies were obtained directly from Es through
an extension of the MUSIC paradigm. This solution, referred
to as Space-Time MUSIC (ST-MUSIC), bypasses traditional
filtering, SCM building and EVD stages. In the sequel, the
STCM signal model at frequency ν is generalized to the non
spherical noise case and the ML ST-MUSIC subspace estimate
is derived from the solution of an inversion problem from Es .
The space-time array response to a pure unit amplitude,
complex sinusoid s (n) = ejνn , impinging from the generic
angle θ, can be written as
xST (n) = aST (ν, θ) ejνn

(16)

where aST (ν, θ), referred to as the space time steering vector
(STSV), has size M P × 1 and is defined by the Kronecker
product
aST (ν, θ) = a (ν, θ) ⊗ eP (ν) = P 1/2 E (ν) a (ν, θ)

(17)

where a (ν, θ) is the classical M × 1 narrow-band array steering vector at frequency ν [20] and eP (ν) =

T
is a vector of size P ×1 [27],
1 e−jν · · · e−j(P −1)ν
[31]. The last equality in (17) defines the M P ×M orthogonal
matrix
E (ν) = P −1/2 IM ⊗ eP (ν)
(18)
as a basis for the subspace spanned by any aST (ν, θ).
To put into evidence the relationship between aST (ν, θ) and
H (θ), let us define the unitary matrix of size P + Ld − 1
h
i
D(ν) = (P + Ld − 1)−1/2 eP +Ld −1 (ν) D⊥ (ν)
(19)
where D⊥ (ν) is the orthogonal complement [36] to
eP +Ld −1 (ν).
After posing sF,d (n, ν) = DH
Ld +P −1 (ν) sd (n) and
H
 (ν, θd ) = H (θd ) DLd +P −1 (ν) =

aST (ν, θd ) H (ν, θd ) D⊥ (ν)

(20)

the response (6) of the d-th source signal after noise whitening
can be rewritten by (16) as
−1
−1
v xST (n) = Rv H (ν, θd ) sF,d (n,ν) =
 R1/2
P E (ν) a (ν, θd ) H (θd ) D⊥ (ν) sF,d (n, ν) .
(21)
To describe the effects of noise whitening on the signal
model at frequency ν, let us define the following full QR
decomposition (QRD)



 R−1
v (ν)
A
(ν)
A
(ν)
P 1/2 R−1
E
(ν)
=
(22)
⊥
v
0

R−1
v

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

where A (ν) is the M dimensional frequency subspace basis
at ν, A⊥ (ν) is its orthogonal complement and R−1
v (ν) is a
square upper triangular matrix of size M , which defines the
following noise whitened narrow-band steering vector of size
M ×1
b (ν, θ) = R−1
(23)
v (ν) a (ν, θ) .
Combining (21), (22) and (23) leads to the decomposition


AH (ν)
R−1
H
v H (ν, θd ) sF,d (ν, n) =
 A⊥ (ν)

(24)
b (ν, θd ) B12 (ν, θd )
sF,d (ν, n) .
0
B22 (ν, θd )
The first M rows contain the sum of the desired, rank one
signal component at frequency ν, characterized by b (ν, θd ),
and of a spectral leakage component transmitted by the
M × (P + Ld − 2) matrix B12 (ν, θd ). The latter term arises
from the finite array aperture and represents the multiple
rank response of a spatially spread, uncalibrated ghost source
[12]. Its strength depends on the spectrum of sd (n) and its
spatial spread on the change of a (ν, θ) with frequency, which
increases with the fractional bandwidth and the array aperture.
Finally, the isolated spectral leakage components can be
observed in the subspace A⊥ (ν) through the transfer matrix
B22 (ν, θd ).
Eq. (24) is immediately generalized to the entire Es (15) as
 H




A (ν)
B11 (ν) B12 (ν) CF 1 (ν)
Es =
(25)
0
B22 (ν) CF 2 (ν)
AH
⊥ (ν)


where B11 (ν) = b (ν, θ1 ) · · · b (ν, θD ) is the M ×
D (ν) narrow-band array transfer matrix at frequency ν of
the D (ν) ≤ D sources with non-zero spectra at the same ν,
combined by the matrix CF 1 (ν) of size D (ν) × η.
the mixing matrix CF 2 (ν) of size
PFinally
D
(P
+
L
× η models the spectral leakage
d − 2)
d=1
observable through the transfer matrices


B12 (ν) =  B12 (ν, θ1 ) · · · B12 (ν, θD ) 
(26)
B22 (ν) = B22 (ν, θ1 ) · · · B22 (ν, θD )
that mark the departure of the convolutional signal model (25)
from the ideal narrow-band one [20].
A. ST-MUSIC as an Inversion Problem
With reference to (25), a consistent narrow-band signal
subspace estimator must asymptotically recover for N/M P →
∞ an M × κ (ν) linearly independent basis Bs (ν) (with
κ (ν) ≤ D (ν) < M ) within the span of B11 (ν) and a full
rank, square subspace weighting matrix Cs (ν) of size κ (ν),
chosen to optimize the statistical accuracy, satisfying the WSF
equation
Bs (ν) Cs (ν) = B11 (ν) C11 (ν)
(27)
for a proper full rank mixing matrix C11 (ν) of size D (ν) ×
κ (ν) [6].
In this formulation κ (ν) is the number of uncorrelated
source signals that are active at ν [42]. As for the SCM case,
spatially coherent (multipath) sources with Ld + Ls,d < P
have a single rank signature in Bs (ν) and the DOA identifiability is subjected to the same limitations [6]. In particular,

5

the steering vectors of the fully coherent arrivals must be
calibrated and Bs (ν) must be analyzed by a multi-dimensional
(e.g., WSF type) DOA estimator [6], [25] or subjected to
spatial or frequency smoothing [8], [12], [43].
As usual, for unambiguous DOA identification, it is required
that D (ν) < M and that any subset of D (ν) steering vectors
b (ν, θ) is linearly independent, at least in a neighborhood of
the source DOAs [20].
A viable method for consistently recovering Bs (ν) from
(25) is to apply a η × κ (ν) inversion weight matrix W (ν) to
the right side of Es to find a M dimensional basis entirely
lying within the span of A (ν):
Es W (ν) = A (ν) Bs (ν) .

(28)

In particular, W (ν) must cancel the out of band Es
components in the subspace A⊥ (ν)
B22 (ν) CF 2 (ν) W (ν) = 0

(29)

to cancel also the spectral leakage in the bin subspace A (ν)
B12 (ν) CF 2 (ν) W (ν) = 0 .

(30)

All involved matrices can be rank deficient under certain
conditions. In general, it is immediate to show that the signal
subspace can be exactly recovered iff
2

|B12 (ν) CF 2 (ν) w|2
2

|B22 (ν) CF 2 (ν) w|2

<∞

(31)

for any non-zero vector w which does not lie in the
intersection of the null-spaces of B12 (ν) CF 2 (ν) and
B22 (ν) CF 2 (ν)4 .
Condition (31) just moves to the frequency domain the
requirement of non-perfect coherency among narrow-band
sources for the applicability of the spatial-only MUSIC [20]
(i.e., spectral vs. spatial coherency), but cannot hamper the
identifiability of spatially coherent sources from Bs (ν).
The structure of (15) reveals that W can cancel out the
spectral leakage from Es separately for each uncorrelated
signal, so Bs (ν) contains linear combinations of all the
steering vectors of active sources. Therefore, the limit STMUSIC equation [22], [31],
EH
v A (ν) Bs (ν) = 0 ,

(32)

herein obtained by projecting both sides of (28) onto Ev is
a sufficient condition for the consistency of a signal subspace
estimate.
DOA estimators starting from approximations of (32) were
considered in the past [4], [13], [42]. However, they employed
non consistent, SCM based noise subspace formulations, basically unable to exploit the source power information and
fruitfully deal with spatially coherent scenarios.
On the contrary, with consistent estimates of Es and Ev
for N/M P → ∞, the sample ST-MUSIC subspace B̂s (ν)
asymptotically matches the ideal narrow-band signal subspace
4 Intersection of these null spaces may not be empty for sources with spectral
nulls within the array bandwidth, but in this case there is not any leakage to
cancel. Infinite (31) indicates that a component does exist at frequency ν and
not elsewhere, so it must be included in B11 (ν).

6

SUBMITTED

model [20] and retains the full flexibility of the WSF approach
in narrowband and broadband scenarios [5], [6], [12] .
A different problem of ST-MUSIC is the unduly attenuation of some sources of interest by (28), evaluated as
B11 (ν) CF 1 (ν) W (ν). In particular, harmonic sources made
up by less than L + P − 1 sinusoids with all frequencies
different from ν or strongly cyclo-stationary sources with cycle
frequency α < N −1 [30] do not satisfy (32) and might be
utterly suppressed. In fact, these sources have a rank deficient
covariance block Sdd in (11).
A basic ST-MUSIC subspace estimator was presented in
[22], starting from a rank-reducing approach applied to Es ,
conceptually similar to TOPS [13]. An approximate STMUSIC for low SNR applications appeared in [31], based
H
on the SVD of the weighted A(ν) Es . These estimators
used biased STCM estimates in finite samples and could not
estimate κ (ν) in a statistically sound way.
For these reasons, in the sequel we set up a ML inversion

problem of Bs (ν), based on the first order (i.e., O N −1/2 )
perturbation model [6], [12], [36], [44] of the sample counterpart of (13)
−H
R−1
= ÊΛ̂ÊH
(33)
v R̂ST Rv
where
R̂ST =

N
X
1
xST (n) xH
ST (n)
N −P +1

(34)

n=P

is the unbiased STCM estimate built from N array output
samples (1). The sample eigenvalues of (33), λ̂k = Λ̂ (k, k)
are ordered for k = 1, 2, . . . , M P in a non-increasing manner.
The corresponding orthonormal sample eigenvectors Ê (:, k)
are partitioned into the sample wide-band signal subspace
Ês = Ê (:, 1 : η), of dimension 0 ≤ η < M P , related to the
dominant λ̂k , and the complementary sample wide-band noise
subspace Êv = Ê (:, η + 1 : M P ), related to the smallest
M P − η eigenvalues, clustered around λv [22], [31].
The EVD version of (33) is preferable since the statistical
identification of a spherical Êv ensures the feasibility of
consistent DOA estimation and reduces the parametrization
and the Mean Square Error (MSE) w.r.t. the true STCM [32].
Under the hypotheses made in Sect. III, two distinct STSs
(2) xST (n1 ) and xST (n2 ) cannot be considered as statistically independent if |n1 − n2 | ≤ P . However, using
independent STSs would require N  M P 2 to get a stable
R̂ST .
In contrast, theoretical arguments [45] and past experience
[22], [27], [31] show that (34) converges only for N  M P
to its asymptotic performance, very similar (except for leakage
issues) to the one of a set of SCM estimates, using the same
N and a DFT length Pw = P , despite the much larger number
of degrees of freedom of the STCM. The following analysis
sheds light on this observation.
V. ML ST-MUSIC S UBSPACE E STIMATOR
The present STCM analysis extends the classical one developed for the SCM [6], [36], [44], [46], points out the
differences and is validated by the simulations in Sect. VII. In
particular, it shows that the ML ST-MUSIC subspace departs

from the sample version of (32) [22] and leads to a generalized
eigen-problem.
A. Sample STCM eigenvectors
The sample bases Ês and Êv are rotated versions of the
true ones in (13) [36] and converge to these for N/M P → ∞
[44]. The (first-order) partitioned rotation matrix is defined by



 
 I + Gss
Gns
(35)
Ês Êv ≈ Es Ev
Gsn
I + Gnn
where Gns = −GH
= −GH
=
sn , Gss
ss and Gnn
H
−Gnn are random perturbation derivatives with L2 norm of
O |RST |2 N −1/2 [44].
The analysis is simplified by recognizing
that blocks

I + Gss and I + Gnn are O N −1/2 approximations to unitary matrices [47]. The invariant subspace Ev is intrinsically
defined up to an arbitrary rotation. The rotation within Ês
was considered in random matrix theory [48] and covariance
shrinkage [32] and is strong for close signal eigenvalues.
While the influence of this rotation on the accuracy bounds
of DOA estimates is unknown, (32) shows that full rank
linear mixing of Es is unessential for estimator derivation,
since a backward transformation is induced on W (ν) by the
numerical optimization on the particular realization of Ês . In
addition, attempts to take into account Gss may face undesired
problems [32].
After these observations, effects caused by Gss and Gnn
can be mostly neglected. Instead the random entries of Gsn
are the most relevant for DOA estimation. They are the finite
sample cosines between Ês and the true Ev , calculated as [6],
[36], [44]
ρ̂sn (p, q)
(36)
Gsn (p, q) =
λq − λv
for p = 1, 2, . . . , M P −η and q = 1, 2, . . . , η. In this equation,
N
P
1
ρ̂sn (p, q) = N −P
yη+p (n) yq∗ (n) is the finite sample
+1
n=P

correlation between the uncorrelated signals yη+p (n) and
H
yq (n), where the generic yk (n) = E(:, k) R−1
v xST (n) is
the signal extracted by the transformed true k-th eigenvector R−H
v E (:, k), considered as a DS beamformer [37] and
referred to as eigenfilter.

Therefore, E {ρ̂sn (p, q)} = 0 + O N −1 , while the covariances among Gsn entries are affected by the temporal
correlation between yη+p (n) and yq (n) [49]. For stationary
signal and noise, the following covariances are computed after
[33], [44]
E {ρ̂sn (p, q) ρ̂∗sn (k, l)} =

N −P +1−τ
E yη+p (n) yq∗ (n) ×
(N −P +1)2
τ =−(N −P )
o
∗
(n + τ ) yl (n + τ )
yη+k
NP
−P

NP
−P

E {ρ̂sn (p, q) ρ̂sn (k, l)} =
N −P +1−|τ | 
E yη+p (n) yq∗ (n) ×
(N −P +1)2

τ =−(N −P )

yη+l (n + τ ) yl∗ (n + τ )} .

(37)

(38)

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

7

Finally, if signals and noise are realizations of zero mean,
circularly complex processes for every τ with zero third-order
moments, it is possible to approximate (37) and (38) as
E {ρ̂sn (p, q) ρ̂∗sn (k, l)} ≈ (N − P + 1)

−1

δpk δql λv λq (39)

E {ρ̂sn (p, q) ρ̂sn (k, l)} = 0

(40)

independently of the exact distributions of signals and noise
[44]. For N  M P , the entries of Gsn approach a Gaussian
distribution by a Central Limit Theorem argument [6]. Then
asymptotically (36) can be approximated as
1/2

Gsn ≈ G0 Γ0

(41)

where G0 is a (M P − η) × η random matrix with i.i.d., zero
−1
mean circular Gaussian entries of variance (N − P + 1)
and
nh
io
1/2
1/2
1/2
Γ0 = diag
(42)
γ0 (1) · · · γ0 (η)
−2

where γ0 (k) = λk λv (λk − λv ) , exactly like in the narrowband case using N − P + 1 independent snapshots [6].
Because of the random rotations Gss and Gnn , G0 is
not observable, but filtered versions of it enter the analysis
as large random Wigner covariance matrices (e.g., GH
0 G0 ).
Wigner matrices converge to a common limit distribution
with Gaussian off-diagonal entries and a semicircle eigenvalue
distribution, regardless reasonable violations of (39) and (40)
[51]. This observation well explains the empirical results of
[22], [27], [31].
Curiously, the largest deviations are expected for very small
γ0 (k) ≈ λv /λk , i.e., at high SNR, and in the presence of
strong temporal correlation of source signals.
Fig. 2. Top: Sample noise STCM eigenvector beampattern versus wavenumber
and angular frequency, normalized to the array bandwidth center for four
pass-band, uncorrelated sources, impinging on the ULA used for simulations
in Sect. VII-B for SNR= 20 dB. Bottom: Average beampattern computed
from the frequency-interpolated least dominant SCM eigenvectors of 64 DFT
bins from the same realization, showing null depth reduction and bias due to
leakage artifacts.

With a large computational effort, some statistics of (37)
and (38) might be estimated from the sample yq (n). However,
a great simplification can be obtained by examining the
properties of STCM eigenfilters.
In fact, the beampatterns of signal eigenfilters in the frequency/DOA space exhibit main-lobes pointing at the regions
of activity of sources and extract convolutive mixtures yq (n)
(q = 1, 2, . . . , η) of source signals plus noise.
On the contrary, as shown in Fig. 2, which realizes Fig.
1 in a computer simulation, noise eigenfilters generally have
beampatterns with very deep broad-band nulls steered toward
the source DOAs [4]. Therefore, they extract a set of noise
dominated signals yη+p (n) (p = 1, 2, . . . , M P − η) that
are almost white and mutually uncorrelated and additionally
uncorrelated with any source dominated signal yq (n) (q =
1, 2, . . . , η), because (33) removes in a Least Squares (LS)
sense [50] the temporal correlation from yη+p (n) for any
τ < P.

B. Normalized signal eigenvectors
Perturbed signal
 eigenvectors in (35) are orthonormal only
up to O N −1/2 terms [44]. From (42), these eigenvectors
exhibit a non-physical singularity for λk → λv , when sample
unit norm eigenvectors only slip in the estimated Êv . To circumvent this problem, in the sequel we will normalize the k-th
perturbed p
signal eigenvector (35) with respect to its expected
L2 norm 1 + cγ0 (k), where c = (M P − η)/(N − P + 1)
is the ratio between the dimension of Ev and the number
of STSs used in (34) and resembles a parameter defined in
random matrix theory[48], [52].
This fixed O N −1 scaling cannot modify the asymptotic
subspace performance for N → ∞ [44]. Even if perturbations
do not likely follow the first order modeling for λk → λv , this
scaling brings out correction terms that taper off the influence
of signal eigenvalues down to zero in the same limit and
partially move marginal eigenvectors into Êv . However, the
corrections are not negligible for typical values of N and Γ,
therefore increasing the robustness to estimation errors of η.
The scaled perturbation model (35) of Ês , neglecting the
inessential rotation induced by Gss , becomes
Ês ≈ Es Γ1/2
+ Ev G0 Γ1/2
s
1/2

(43)

where Γs and Γ1/2 are the diagonal square roots [36] of
−1
−1
Γs = [Iη + cΓ0 ]
and Γ = Γ0 [Iη + cΓ0 ] . In addition,
0 < cΓ (k, k) < 1 and (41) can be recovered for c → 0.

8

SUBMITTED

C. Signal subspace rank estimation for the STCM
From the previous arguments, it is clear that the rank η of Es
is not directly related to the number of uncorrelated impinging
signals, as in the narrow-band case [23], but it represents
the number of STCM signal components that exceed the
noise level λv and carry significant information about source
parameters. In essence, the problem is to find a spherical
noise subspace centered on a consistent O N −1/2 estimate
of λv . This goal involves the STCM eigenvalue statistics, still
affected by the snapshot dependency issue [49].
Under the same assumptions made in Sect. V-A, the perturbative model for the sample STCM eigenvalues λ̂k starts
N
−1 P
2
from the approximation λ̂k ≈ (N − P + 1)
|yk (n)|
n=P
n o

[36], [44]. It is deduced that E λ̂k = λk + O N −1
n o

for k = 1, 2, . . . , η and E λ̂k = λv + O N −1 for
k = η + 1, η + 2, . . . , M P .
For large N , the covariance between λ̂k and λ̂l , corresponding to distinct λk 6= λl , as well as the variance of isolated
signal λ̂k , approach [44], [49]


NP
−P
N −P +1−|τ |
×
cov λ̂k , λ̂l '
(N −P +1)2
τ
=−(N
−P
)
(44)
n
o
2
2
E |yk (n)| |yl (n + τ )| − λk λl .
In particular, pairs of signal and noise sample eigenvalues
are almost uncorrelated. However, (44) indicates that the
variance of signal STCM eigenvalues is much higher in the
dependent snapshot case, but still generally lower than in the
SCM case.
As regards the sample noise eigenvalue statistics, the consistency of the classical noise variance estimate [23]
λ̂v =

MX
P −η
1
λ̂k
MP − η

(45)

k=η+1

follows from the convergence of (34) to the true STCM under
the given assumptions on the choice of P [45].
In addition, the robust censoring estimator of η proposed
ensemble variance
n in [12] showed that the sample
o
var λ̂η+p ; p = 1, 2, . . . , M P − η
of (33) is always close
to cλ2v for Gaussian noise, as in the independent snapshot
case [52]. This result confirms that the extracted noise signals
yη+p (n) behave as white and mutually uncorrelated for any
τ of interest, in agreement with the analysis of [45] for
more general noise distributions. However, a small increase
of this ensemble variance was observed at low SNR and was
interpreted as a symptom of increasing leakage of temporally
correlated signals into Êv .
Therefore, by Chebyschev’s inequality, the lower threshold
√
for signal eigenvalues should be a small multiple of cλ̂v
above λ̂v to avoid missing the small components originated
from the fading tails of the impulse responses hm (θd ) [22].
It follows that any optimal value of η shrinks with the SNR
and this phenomenon may impact on the DOA identifiability
by (32), since a reduced dimension W (ν) may not consistently recover the full span of Bs (ν) anymore.

Information Theoretic criteria [23], [52] developed for the
SCM are questionable for the STCM, because signal and noise
processes exhibit very different and a priori unknown numbers
of effective observations by (44) and the likelihood form is
different for each signal distribution.
However, the Gaussian log likelihood ratio [23] for independent observations essentially relies on the noise eigenvalue
spread, that changes little under the assumption of Sect. V-A,
even for many non-Gaussian noise distributions [45]. These
arguments support the provisional applicability of unmodified
existing SCM rank estimation criteria to the STCM. In a
refined model, the dependent observation issues might be
handled by assuming a smaller N̄ < N − P + 1 number
of observations [49] in the relevant noise subspace, but a
quantitative study is beyond the scope of this work.
In particular, the AIC [23], which takes into account only
the validation risk [35] of adding a new signal eigenvector,
performed better, especially at low SNR, in comparison with
the non parametric approach of [12], the BIC [23] and the
modified AIC of [52], that evidenced clear mismatches of their
additional parameters. In particular, despite of consistency
claims and stable η estimation at high SNR, the BIC exhibited
a disappointing detection threshold at low SNR, about 10 dB
higher than the AIC one. Therefore, only the AIC will be
included in the simulations of Sect. VII.
D. ML ST-MUSIC inversion
The main issue with earlier versions of ST-MUSIC [22],
[31], based on the SVD of ÊH
s A (ν), was the statistical
instability of the signal subspace weighting, which generalized
existing weighted MUSIC schemes [25], [33]. In particular, the
empirical estimate Γ̂ of the large Γ led to signal cancellation
at high SNR, as in adaptive beamforming [28]. Heuristic
regularization of Γ̂, based on [32], stabilized the subspace
estimate at the expense of the sample DOA bias, partially
spoiling the goal of consistency.
A true ML formulation circumvented these issues by rewriting the sample version of (28) at frequency ν, after inserting
(43), as


Ês W (ν) = A (ν) B (ν) + Ev G0 Γ1/2 W (ν) + o N −1/2
(46)
where B (ν) is the unknown candidate signal subspace basis
at (ν) for Bs (ν) in (27).
It is assumed that η has been estimated as described in
Sect.V-C. The true Γ and W (ν) are assumed initially known.
In the sequel, for conciseness we will set NP = N − P + 1
and drop the reference to ν, e.g., A = A (ν), B = B (ν) and
W = W (ν).
The sample equation error Ev G0 Γ1/2 W of (46) lies in
the true Ev . One of the instability sources was found in the
classical replacement of Ev by the sample Êv , that called for
a more sophisticated route. Projecting (46) onto Êv and Ês
leads to


1/2
−1/2
ÊH
W = ÊH
(47)
v Ev G0 Γ
v AB + o N




1/2
−1/2
Iη − Γ1/2 GH
W = ÊH
(48)
0 G0 Γ
s AB + o N

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

w.r.t. W, leading to the asymptotically unbiased estimate

and


1/2
1/2
ÊH
W = Γ1/2 GH
W + o N −1/2
s Ev G0 Γ
0 G0 Γ



.
(49)
Equations (47) and (49) together characterize a row rotation
of Ev G0 Γ1/2 W, while (48) can be used to estimate W. After
(41), G0 Γ1/2 W can be modeled without loss of generality as
a stack of M P − η independent, zero mean, Gaussian circular
observations of dimension M with covariance NP−1 WH ΓW ,
characterized by the following negative log likelihood for
given W and Γ

(53)


 −1
where Φ = Iη − 2cΓ + c2 + cηNP−1 Γ2
(Iη − cΓ) is a
real-valued, diagonal matrix of size η.
The candidate solution for B can be found from (52)
through the Generalized SVD (GSVD) [36], which leads in
a square root fashion to the decompositions
H

2 H
2 H
2
Πε = AH Êv ÊH
(54)
v A+c A Ês ΦΓ ΦÊs A = FΣε F
H

(50)

(0)

Πε = BH AH Êv ÊH
v AB+
H 1/2 H
1/2
W Γ G0 G0 ΓGH
W
0 G0 Γ

H

Ŵ = ΦÊs AB

Πs = AH Ês ΦΓΦÊs A = FΣ2s FH

ϕ (M ) = M (M P − η) ln (π) +

(M P − η) lnh det NP−1 WH ΓWi +

−1
(0)
NP trace Πε WH ΓW
where

9

.

(0)

The last term of Πε can be replaced by its expected value
under the Gaussian i.i.d. approximation on G0 , leading to

(55)

where Σ2s and Σ2ε are diagonal square matrices of size M with
non-negative diagonal entries, linked by Σ2ε + Σ2s = IM , and
F is a complex valued square matrix of size M . By posing
B = F−H [36], (52) can be simplified as

ϕ (M ) = M (M P − η) ln πNP−1 +
(56)
(M P − η) ln det Σ2s + cM trace (Γ) +
.
NP trace Σ2ε Σ−2
s

Not every column of B is acceptable as a basis vector for the
(0)
H H
H
AB+
=
B
A
Ê
Ê
Π
signal
subspace, since ϕ (M ) has small contributes only along
v v
ε

(51)
the κ columns of B characterized
WH c2 Γ2 + cΓ · NP−1 trace (Γ) W .
 2 by very small generalized
2
eigenvalues
µ
=
Σ
(k,
k)
Σs (k, k) ' c.
k
−1
ε
The term WH ΓW
in (50) can be interpreted as the
In fact, plugging (43) into (54) and (55) and taking the
high resolution Capon estimate [19] of the power spectral expected value over G0 , leads to the following O N −1
density (PSD) at frequency ν of a set of signals characterized approximations
by the covariance Γ−1 and the steering vector matrix W. As a
o
n


consequence, signals components at frequencies different from
E {Πε } ' IM − E AH Ês Iη − c2 ΦΓ2 Φ ÊH
s A
(57)
ν, such as spectral leakage products, are suppressed as in MV
' AH Es cΓ (Iη + cΓ)
EH
s A+

−1
2 3 2
H
H
beamforming [28]. The same term in (51) compensates for the
1 − NP trace Γ − c Γ Φ A Ev Ev A
inclusion of marginal
signal eigenvalues close to λ̂v .

cΓ) EH
E {Πs } ' AH Es Γ (I
The O N −2 term c2 WH Γ2 W in (51) cannot influence
s A+
η +
(58)
−1
2 2
H
NP trace Γ Φ A Ev EH
the asymptotic performance, but it is the key for regularizing
v A
(0)
(50) at any SNR by filling the numerical rank of Πε and where some non-negligible terms from the Taylor series excompensating for any rotation of the sample W. Its exact pansions of (53), (54) and (55) w.r.t. c are retained for clarity
scaling mildly depends on the assumed distribution of GH
0 G0 ,
and future use.
but the given calculus for the Gaussian case is sufficient for
Thus, on the average, for N → ∞ the generalized eigenregularization purposes.
equation corresponding to (54) and (55) [36]
Plugging (51) into (50) and simplifying leads to

E {Πε } B (:, k) ' µk E {Πs } B (:, k)
(59)
ϕ (M ) = M (M P − η) ln πNP−1 +
H
(M
hP − η) ln det W ΓW + cM trace
 (Γ) + NP× i is asymptotically satisfied by an eigenvalue µk = c of multi−1
H H
H
2
H 2
trace B A Êv Êv AB + c W Γ W WH ΓW
. plicity κ. The corresponding κ eigenvectors, collected in the
H
H
(52) M × κ matrix Bs , asymptotically satisfy A Ev Ev ABs ' 0
There is some freedom in choosing an estimate Ŵ of W and therefore (32).
In addition, the null-space of E {Πε − cΠs } is independent
from (48). Early weighted forms of MUSIC DOA estimators
of
N for sufficiently large samples, demonstrating that the
[25], [33] obtained null spectra resembling the last term of
span
of the sample Bs (i.e., B̂s ) provides an asymptotically
H
(50) for the implicit choice Ŵ = Ês AB. However, errors
unbiased
in Ŵ may destroy the Capon spectral estimate of (52) [28].
 ST-MUSIC signal subspace estimate up to at least
−1
O
N
terms.
Therefore, we sought for linear estimates of W, independent
This
claim
is not evidently shared by earlier
ST-MUSIC
of G0 , that minimize both the residual error of (48) and the

−1
estimators
[22],
[31]
that
retain
O
N
signal
subspace
variance of Ŵ, according to (43). In particular, the form of
n
o
H
H
the random term of (43) suggested to minimize the expected components in E A Êv Êv A , related to the source specvalue over G0 of the weighted LS functional
tra. In fact, the different asymptotic bias at various frequencies
h

 i2
is another source of DOA estimate instability at high SNR in
1/2 H
Γ−1/2 ÊH
G0 G0 Γ1/2 W
wide-band applications without a regularized Γ̂.
s AB − Iη − Γ
F

10

SUBMITTED

However, beside these finite sample effects, the asymptotic
variances of the various ST-MUSIC subspaces are identical, as
for the narrow-band MUSIC DOA estimator [25], [33], [48].
For the leakage components to be discarded, instead, Σ2ε (k, k) is generally close to one according
to
h (32), so µk is comparablei to the Capon PSD
H

H

−1

B(:, k) AH Ês ΦΓΦÊs AB (:, k)
, which is generally
much greater than c [52]. However, in rare scenarios, leakage
residuals might generate some µk of order c, interpreted as
extremely weak ghost sources, that have a vanishing impact
for N → ∞ and can be therefore harmlessly included in B̂s .
Finally, the M − κ columns of B, corresponding to the
discarded µk  c, are collected in the matrix B̂ε .

The estimates of κ abruptly shrink at very low SNR,
as in the sample SCM case [52]. However, the eigenvectors discarded by (61) and (62) are certainly dominated by
whitened leakage residuals, are numerically wobbly and must
be excluded from B̂s .
On the other hand, (60) involves a small ratio between the
number of observations M P − η and κ < M and the risk for
DOA estimation of including marginal µk s originated by false
alarms is low. Therefore, the consistency claims of BIC for
(M P − η)/κ → ∞ are weak and the slightly more permissive
AIC (61) may alleviate the κ shrinkage issue at low SNR and
was adopted in simulations.
F. Optimal ST-MUSIC Subspace Weighting

E. ST-MUSIC Rank Selection by AIC and BIC
Due to the high spectral variability of wide-band sources,
the signal subspace rank κ generally changes with frequency
and should be estimated from data. The heuristic rules based
H
on the singular value magnitude of A(ν) Ês , used in earlier
ST-MUSIC approaches [22], caused ambiguities.
On the contrary, the proposed ML formulation admits the
use of Information Theoretic Criteria [23], [34], [35] for
estimating κ. To this purpose, sample µk s are ordered in a non
decreasing manner5 and a set of nested models with reduced
rank K is compared for 0 ≤ K < M . The log-likelihood of
each model is written as

ϕ (K) = K (M P − η) ln πNP−1 +
(M P − η) lndet Σ2s (1 : K, 1 : K) + cKtrace (Γ)
 +
NP trace Σ2v (1 : K, 1 : K) Σ−2
s (1 : K, 1 : K)
(60)
where, by a continuity argument, ϕ (0) = 0. The estimates of
κ by AIC and BIC are readily [23], [34], [35] obtained as

Optimal DOA estimation in the WSF framework [6],
[25] entails the calculus at each frequency ofn interesto of
−
the spatial perturbation covariance Ξε = E B̂s B̂H
s
n o n
o
n o
H
E B̂s E B̂s , where E B̂s = Bs from (59) and of the
optimal subspace weighting matrix Cs nof size
o κ × κ, which
whitens the perturbation of Bs onto E B̂ε = Bε .

This task requires the O N −1/2 perturbative analysis of
(59), adapted and simplified from [36], [44], [53]. The eigenvector perturbation Ḃs is herein decomposed for convenience
as
Ḃs = B̂s − Bs = Bs Ys + Bε Yε
where Iκ + Ys describes a random rotation of Bs , since all
signal eigenvalues µk → c asymptotically.
Differentiating (59), Yε asymptotically leads to


− Π̇ε − cΠ̇s Bs ' E {Πε − cΠs } Bε Yε
(63)
where, after tedious calculus using (43), (53), (54), (55)

H
Π̇ε ' Πε − E {Πε } ' AH Es [−Iη + 0.5cΓ] Γ1/2 GH
0 Ev A
H
1/2
H
+A Ev G0 Γ [−Iη + 0.5cΓ] Es A
0≤K≤M −1
(64)
H
κ̂BIC = arg min {2ϕ (K) + K (2M − K + 1) ln (M P − η)}
A
Π̇s ' Πs − E {Πs } ' AH Es Γ3/2 GH
E
v
0
(65)
0≤K≤M −1
+AH Ev G0 Γ3/2 EH
s A.
(62)
since the number of free real parameters for the model order
Inserting (57), (58), (64) and (65) into (63) and taking into
K is M 2 for the preliminary whitening GSVD transformation account from (59) that EH
v ABs = 0, we get
(i.e., an unessential constant overhead), plus K (2M − K) for
H
AH Ev G0 Γ1/2 (Iη + 0.5cΓ)
s '

 HEs AB
the rank K subspace B̂s [36], plus K for the generalized
(66)
−1
2
H
1
−
N
trace
Γ
+
cΦΓ
Φ
A
E
E
AB
6
v
ε Yε .
v
P
eigenvalues .
Either (61) and (62) implicitly define an environment deTaking expectations on both sides yields E {Yε } = 0,
pendent upper threshold Tµ > c for a valid µk . Their detailed because E {G0 } = 0, confirming the absence of asymptotic
analysis is complicated by the statistical interaction with the bias. In addition, right multiplying both sides by the unknown
prior estimation of η and is left out of the scope of this paper. Cs , neglecting terms independent of ν, using (58) and solving
However, since c  1 and µk  c for leakage components, for Bε Yε Cs asymptotically yields
(61) and (62) furnished reliable and almost indistinguishable

†
H
H
Ξε ∝ Bs E Yε Cs CH
= AH Ev E H
s Yε  Bs
v A ×
detection results in simulation for both κ and DOA estima2
NP−1 trace Cs CH
.
tion. This result further supports the findings and the basic
s E Σs (1 : κ, 1 : κ)
(67)
assumptions made in the STCM analysis.
 2
−1
Since from (59) E
Σ
(1
:
κ,
1
:
κ)
'
(1
+
c)
I
,
minκ
s

 2


5 The quantity ln Σ2 (k, k) + c−1 Σ2 (k, k) Σ2 (k, k) increases monotonimization of trace Cs CH
s E Σs (1 : κ, 1 : κ) , subject to
s
ε
s
det Cs CH
= 1, asymptotically yields [6]
ically with µk and Σ2ε (k, k).
s

κ̂AIC = arg min {2ϕ (K) + 2K (2M − K + 1)}

(61)

6 The columns of B are orthonormal after the whitening and µ , Σ2 (k, k)
k
ε
and Σ2s (k, k) are all linked one-to-one.

Cs = Iκ

(68)

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

or any other unitary
(59) a consistent estimate
† matrix.−2From
of AH Ev EH
is B̂ε Σκ B̂H
v A
ε , where
Σ2κ = (1 + c) Σ2ε (κ + 1 : M, κ + 1 : M ) − cIM −κ .
Such an error covariance matrix confined within the B̂ε
subspace of dimension M − κ can cause numerical troubles
to WSF estimators especially during the coarse DOA initialization. So the sample Ξ̂ε can be completed by adding
a rather arbitrary Hermitian matrix spanning the orthogonal
complement Qε of B̂ε . In particular, the choice


H
H
Ξ̂ε = κ B̂ε Σ−2
(69)
κ B̂ε + ζε Qε Qε


−1
H
where ζε = (M − κ) trace B̂ε Σ−2
κ B̂ε , effectively minimized the distance Ξ̂ε − κζε IM

F

in most cases with ζε ≈ 1.

However, in the presence of strong leakage Ξ̂ε can still be illconditioned, hampering for instance the statistical efficiency
of common focusing schemes [21].
G. DOA estimation from the ML ST-MUSIC signal subspace
Given the associate Ξ̂ε (ν) and Cs (ν), the sample ML
ST-MUSIC subspace B̂s (ν) at frequency ν can replace the
SCM counterpart in any subspace based narrow-band DOA
estimator. For instance, the optimal WSF estimator minimizes
over the DOA parameter set Θ = {θ1 , . . . , θD } the functional
h
i2
1
b ε (ν)− 2 B (Θ, ν) C − B̂s (ν) Cs (ν)
Θ̂W SF = arg min Ξ
Θ,C

F

(70)
where B (Θ, ν) is the tentative steering matrix for B11 (ν) in
(25), Cs (ν) = Iκ from (68) and C is an unknown full rank
mixing matrix [6].
The extension of (70) to multiple frequencies is straightforward [5]. It allows consistent and asymptotically efficient
DOA estimation in a Gaussian scenario at ν, thanks to the
ML formulation of the ST-MUSIC subspace and of the WSF
minimum DOA variance property [6]. Consistency holds for
any fourth-order bounded scenario in the absence of other
model mismatches [12], [24], [26].
The estimator (70) analyzes spectral slices of the STCM
[19] with a very narrow effective bandwidth around ν, comparable to the width of a null of (50), virtually extrapolating the
STCM and eliminating the effects of steering vector changes
with frequency, beside spectral leakage.
Following this argument, working on P equi-spaced frequencies may not optimally exploit the source information by
ML ST-MUSIC, since some strong narrow-band components
might be missed. In particular, (15) indicates that an optimal
ST-MUSIC would require the coherent analysis of at least
P + Ld − 1 optimally selected points of the full z plane (i.e.,
by using a complex ν) [41], [42].
In addition, since the consistency of signal subspace fitting
does not require the exact
 specification of Ξε (ν) and Cs (ν)
[6], [12], the O N −1/2 estimates (68) and (69) are adequate.
However, obtaining robustness to steering vector mismatches
at ν may require different choices of (68) and (69) [12], [24],
[26].

11

In contrast, with reference to (25), the energy from adjacent frequencies in the span of B11 (ν) (i.e., the component
B11 (ν) B†11 (ν) B12 (ν) CF 2 (ν)) is retained by the SCM and
canceled by the ML ST-MUSIC, which may give to Fourier
methods a detection advantage at low SNR with a coarse
frequency binning.
The overall statistical efficiency of the ML ST-MUSIC
subspace, which, as the SCM, is a reduced statistic w.r.t.
the STCM, remains an open question. However, working in
an ideal narrowband scenario (i.e., after posing A = IM ),
the weightings of the ML ST-MUSIC and of the optimal
narrowband subspaces [6] are proportional.
The STCM allows more flexible strategies in non-stationary
environments. For instance, a STCM can be built from many
disjoint blocks of few STSs each to locate intermittent, elusive
sources without sacrificing the spectral resolution and the
effective SNR as in the binning approach [18].
VI. C OMPUTATIONAL ANALYSIS
The ML ST-MUSIC computing cost is largely dominated
by the building and the EVD of the STCM. In particular, the
STCM computation requires a bulk of M 2 N P multiply and
accumulate operations using the apparatus developed for fast
AR identification [50]. The dominant EVD cost of the sample
STCM is about 6M 3 P 3 complex flops. Viable alternatives
include fast Lanczos algorithms [36].
This effort must be compared to the about 6M 3 Pw + N M 2
flops required by the SCM approach for the same tasks. To
give an idea, in simulations using MATLAB, the full WAVES
using ML ST-MUSIC ran in about 1.0 s, while the SCM based
WAVES in about 25 ms for N = 6400, M = 8, P = Pw = 64
and 33 analyzed frequencies on a PC machine featuring an
Intel Core i7-6700K processor running at 4.2 GHz and 32
GB RAM. However, in most real-time settings, the marginal
cost of STCM processing is comparable with the overhead of
coarse DOA initialization techniques using high performance
wide-band beamformers [15], [16], [54], orthogonal matching
pursuit for WSF [5], [6], or refined focusing schemes [10],
[12], [21], [54].
VII. C OMPUTER SIMULATIONS
The performance of ST-MUSIC was assessed by MonteCarlo time domain simulations7 of a ULA with M = 8 omnidirectional sensors, equi-spaced by 0.5 wavelengths at the
center frequency ν0 . The sensor bandwidth ranged from 0.6ν0
to 1.4ν0 . The background Gaussian noise was temporally and
spatially white (Rvv = IM P ) and the source SNR was referred
to each array element.
For each SNR value, 1000 independent Monte Carlo trials
were run, collecting N = 6400 sensor snapshots sampled at
−1
T = (0.8ν0 )
s and using P = 64 (NP = 6337) for the
STCM and an un-windowed DFT of Pw = P points over
7 It is worth noting that inconsistency effects cannot be observed if wideband array outputs are simulated in the frequency domain by a set of
independent narrow-band signals.

12

Fig. 3. Plot of the DOA sample standard deviation versus the SNR for the
source impinging from 37◦ , observed at the single frequency ν0 , obtained
by MODE [25] applied to different signal subspace estimates, compared with
the narrow-band SCM square root CRB for Ns = 100.

Ns = 100 non overlapping, consecutive data segments for
SCM based approaches8 [8], [12], [25].
Four sources emitting various kinds of signals were placed
at the azimuth angles 8◦ , 13◦ , 33◦ and 37◦ , referred to the
array broadside. Performance comparison employed state of
the art WSF type DOA estimators, namely MODE [25] and
WAVES [12] followed by MODE, in turn fed by sets of signal
subspaces drawn by ML ST-MUSIC and by SCM, optimally
weighted for finite sample errors according to [6], [12],
(68) and (69). Since DOA variance and bias have different
dominant causes (finite sample errors and model mismatches,
respectively [21]), they were separately analyzed.
A. Single frequency analysis
This experiment tested the consistency and the basic accuracy of the subspace estimates. Sources radiated equi-powered
and uncorrelated wide-band Gaussian noise, filtered in the
band (0.725ν0 , 1.275ν0 ) (see Fig. 2). MODE (i.e., a root
version of (70)) was applied in turn to the ML ST-MUSIC
subspace estimate at ν0 , to the corresponding SCM estimate
at bin 0 and to the reference subspace drawn from a classical
narrow-band SCM [20] built with Ns = 100 independent
snapshots, whose Cramer Rao bound (CRB) for DOA was
also calculated. For the STCM, η = 200 was selected from
the median of the BIC estimates at the highest SNR of 40 dB.
The ST-MUSIC subspace rank κ was selected online by AIC
(61). The signal subspace rank of the sample SCM and the
number of sources searched by MODE were instead fixed at
four.
DOA standard deviation and bias are displayed in Figs. 3
and 4 versus the SNR for the source impinging from 37◦ , the
most critical one for CRB and steering vector change with
frequency. For all algorithms, the low SNR threshold was at
8 Setting P
w = P allowed the use the same focusing matrices in Sects.
VII-B and VII-D. Due to the previous discussion, these choices may not be
optimal for each estimator and scenario.

SUBMITTED

Fig. 4. Plot of the DOA sample bias versus the SNR for the source impinging
from 37◦ , observed at the single frequency ν0 , obtained by MODE applied
to different signal subspace estimates.

around 5 dB. The DFT-generated subspace was clearly nonconsistent. The ML ST-MUSIC extracted more information
than the classical SCM subspace for Ns = 100, obtaining
lower DOA variance, coupled with negligible and similar bias.
These results support the consistency and the superior accuracy
potential of STCM and ML ST-MUSIC and can be directly
extended to multi-frequency WSF estimators [5].
B. Wide band focusing
The test ran in the same scenario, but focusing the center
frequencies of the 33 bins in the band (0.8ν0 , 1.2ν0 ) onto
ν0 using WAVES [12], followed by MODE. In this case,
DOA consistency was basically hampered by the focusing
bias [13], [21] and the test emphasized the quality of the
information combination across frequencies and the robustness to model errors. Unitary focusing matrices were used,
assuming that preliminary beamforming estimated two DOA
clusters centered at 10.50◦ and 35◦ with focusing sectors as
in [8], [12]. [8]. However, the number of focusing angles was
increased to twelve (6.7◦ , 8.6◦ , 10.50◦ , 12.3◦ , 14.1◦ , 31◦ ,
32.3◦ , 33.6◦ , 35◦ , 36.3◦ , 37.6◦ and 39◦ ) for better accuracy
of array interpolation. In fact, the unitary Procrustes DOA
interpolation problem [8], [9], [36] over the sector ensemble S
yields
a focusing matrix equal to the orthogonal polar factor of
R
H
b (ν0 , θ) b(ν, θ) dθ. The original angle choice of [8] did
θ∈S

not adequately cover the effective sector rank and approximate
this integral [55].
The focusing accuracy was still unsatisfactory for our
purposes at high SNR and forced us to set the WAVES
rank to four. As an additional measure, a cascaded, nonunitary sector interpolation matrix [56], unique for all bins,
remapped the principal singular vector of each set formed by
the 33 focused steering vectors at the same DOA onto the
corresponding ULA steering vector at ν0 [57]. This heuristic
correction substantially reduced the bias and still satisfied
Fisher optimality criteria [8].

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

Fig. 5. Plot of the DOA sample standard deviation versus the SNR for
the source impinging from 37◦ , obtained by WAVES [12], focused on 33
bin frequencies, followed by MODE applied to different signal subspace
estimates, and compared with the narrow-band SCM square root CRB and
the WAVES PSM bound computed for Ns = 100.

13

Fig. 7. Plot of the DOA sample standard deviation versus the SNR for
the coherent source impinging from 37◦ , obtained by WAVES, focused on
33 bin frequencies, followed by MODE+WSF, applied to different signal
subspace estimates and compared with the narrow-band SCM square root
CRB computed for Ns = 100.

high SNR, further complicated by the ghost source detection
when using AIC, and was also heavily biased at every SNR.
The DOA variance of the WAVES with ML ST-MUSIC subspaces appeared slightly degraded at extreme SNRs w.r.t. the
single frequency case. The most likely causes are the focusing
errors and the Fisher information loss [8] due to the suboptimal
averaging of Ξ̂ε (ν) across the pass-band, performed by the
chosen fixed focusing scheme [12].
C. Wide Band Coherent Sources

Fig. 6. Plot of the DOA sample bias versus the SNR for the source impinging
from 37◦ for the various subspace estimates.

The CRB and the SCM based WAVES bound under the
perfect subspace mapping (PSM) assumption [9], [12] were
computed for Ns = 100 independent snapshots drawn from
each DFT bin. Moreover, beside the cases of fixed selection
for η = 200 and the SCM rank, two further experiments were
included, where these quantities were all estimated online by
AIC [23].
Sample results were displayed in Figs. 5 and 6 for the
source at 37◦ . The standard deviation plot of Fig. 5 revealed
a slight advantage of SCM over ML ST-MUSIC subspaces
at low SNR, with a threshold for reliable detection slightly
higher for AIC driven rank estimates. The difference appears
essentially due to a slightly higher (not removed) outlier rate
of ML ST-MUSIC at extremely low SNR, because the AIC
detected rank κ shrunk at several frequencies in some trials.
Instead, the overall insensitivity of ML ST-MUSIC to the AIC
selection of η was remarkable, as expected.
The SCM based WAVES exhibited the usual breakdown at

The test was repeated with similar results after replacing
the sources at 13◦ and 37◦ with delayed replicas of those at
8◦ and 33◦ with delays 0.45T and 0.35T respectively [22].
In these trials, there was not any obvious choice for η and
the SCM rank, that were selected by AIC. MODE had local
convergence issues and its DOA estimates were refined by
WSF Newton iterations. Results are shown in Figs. 7 and 8.
The PSB was not shown since it is ambiguous for coherent
arrivals, depending on the focusing strategy.
D. Wide Band Autoregressive Sources
The previous test was repeated after replacing the white
sources with independent 16-th order autoregressive (AR)
sources, whose models were drawn from live recordings at an
airport, assuming ν0 = 2π × 103 rad/s. The source correlation
length ranged from 30T to 50T . The SNR was referred to the
driving Gaussian noise variance, equal for all sources. The
STCM rank estimated by BIC at high SNR was η = 199 in
this environment.
Sample results were displayed in Figs. 9 and 10 for the
source at 37◦ and confirmed previous findings.
VIII. C ONCLUSION
A consistent signal subspace estimator at arbitrary frequencies (ML ST-MUSIC) has been developed by exploiting the

14

Fig. 8. Plot of the DOA sample bias versus the SNR for the coherent source
impinging from 37◦ for the various subspace estimates.

SUBMITTED

finite rank representation of wide-band sources in the signal
subspace of the STCM through a ML inverse fitting, based
on the MUSIC paradigm. AIC and BIC estimators of the
number of sources active at each frequency of interest have
been naturally derived in this framework. The ML ST-MUSIC
approach radically circumvents the frequency spread issues of
classical binning and is theoretically supported by the asymptotic nulling of the spectral leakage, identified as the source
of the statistical inconsistency. Optimal weighting for finite
sample errors allows the use of ML ST-MUSIC subspaces
in any narrow- or wide-band subspace fitting DOA estimator.
Simulations support the consistency, the superior accuracy
(evident at high SNR) and the robustness to mismatches
(evident in automatic model selection) of DOA estimation
derived from ML ST-MUSIC subspaces w.r.t. the SCM based
counterparts.
Finally, the ML ST-MUSIC derivation raises the delicate
question if binning approaches are theoretically justified even
in narrow-band scenarios. In fact, since it is always Ld > 1
in physical arrays, a STCM with low order P should be used
in order to guarantee consistency.
R EFERENCES

Fig. 9. Plot of the DOA sample standard deviation versus the SNR for
the AR source impinging from 37◦ , obtained by WAVES, focused on 33
bin frequencies, followed by MODE, applied to different signal subspace
estimates.

Fig. 10. Plot of the DOA sample bias versus the SNR for the AR source
impinging from 37◦ , for the different signal subspace estimates.

[1] P. M. Schultheiss and H. Messer, “Optimal and suboptimal broad-band
source location estimation,” IEEE Trans. on Signal Processing, vol. 41,
no. 9, Sept. 1993.
[2] J. F. Bohme, “Source-parameter estimation by approximate maximum
likelihood and nonlinear regression,” IEEE Journal of Oceanic Engineering, vol. 10, no. 3, pp. 206–212, July 1985.
[3] M. Doron, A. J. Weiss, and H. Messer, “Maximum-likelihood direction
finding of wide-band sources,” IEEE Trans. on Signal Processing,
vol. 41, no. 1, pp. 411–414, Jan. 1993.
[4] K. M. Buckley and L. J. Griffiths, “Broad-band signal-subspace spatialspectrum (BASS-ALE) estimation,” IEEE Transactions on Acoustics,
Speech, and Signal Processing, vol. 36, no. 7, pp. 953–964, Jul 1988.
[5] J. Cadzow, “Multiple source location - The signal subspace approach,”
IEEE Trans. on Acoustics Speech and Signal Processing, vol. 38, no. 7,
pp. 1110–1125, July 1990.
[6] M. Viberg, B. Ottersten, and T. Kailath, “Detection and estimation in
sensor arrays using weighted subspace fitting,” IEEE Trans. on Signal
Processing, vol. 39, no. 11, pp. 2436 –2449, Nov. 1991.
[7] H. Wang and M. Kaveh, “Coherent signal-subspace processing for the
detection and estimation of angles of arrival of multiple wide-band
sources,” IEEE Trans. on Acoustics, Speech and Signal Processing, vol.
ASSP-33, no. 4, pp. 823–831, Aug. 1985.
[8] H. Hung and M. Kaveh, “Focussing matrices for coherent signalsubspace processing,” IEEE Trans. on Acoustics Speech and Signal
Processing, vol. 36, no. 8, pp. 1272–1281, Aug. 1988.
[9] M. Doron and A. Weiss, “On focusing matrices for wide-band array
processing,” IEEE Trans. on Signal Processing, vol. 40, no. 6, pp. 1295–
1302, June 1992.
[10] T.-S. Lee, “Efficient wide-band source localization using beamforming
invariance technique,” IEEE Trans. on Signal Processing, vol. 42, no. 6,
pp. 1376–1387, June 1994.
[11] S. Valaee and P. Kabal, “Wide-band array processing using a two-sided
correlation transformation,” IEEE Trans. Signal Processing, vol. 43,
no. 1, pp. 160–172, Jan. 1995.
[12] E. D. Di Claudio and R. Parisi, “WAVES: Weighted average of signal
subspaces for robust wideband direction finding,” IEEE Trans. on Signal
Processing, vol. 49, no. 10, pp. 2179–2191, Oct. 2001.
[13] Y.-S. Yoon, L. M. Kaplan, and J. H. McClellan, “TOPS: New DOA
estimator for wideband signals,” IEEE Trans. on Signal Processing,
vol. 54, no. 6, pp. 1977–1989, June 2006.
[14] N. Bienati and U. Spagnolini, “Multidimensional wavefront estimation
from differential delays,” IEEE Trans. on Geoscience and Remote
Sensing, vol. 39, no. 3, pp. 655–664, Mar. 2001.
[15] J. Krolik and D. Swingler, “Multiple wide-band source location using
steered covariance matrices,” IEEE Trans. on Acoustics Speech and
Signal Processing, vol. 37, no. 10, pp. 1481–1494, Oct. 1989.

DI CLAUDIO et al.: SPACE TIME MUSIC: CONSISTENT SIGNAL SUBSPACE ESTIMATION FOR WIDE-BAND SENSOR ARRAYS

[16] E. D. Di Claudio and R. Parisi, “Robust ML wide-band beamforming
in reverberant fields,” IEEE Trans. on Signal Processing, vol. 51, no. 2,
pp. 338–349, February 2003.
[17] D. P. Wipf and B. D. Rao, “An empirical bayesian strategy for solving
the simultaneous sparse approximation problem,” IEEE Transactions on
Signal Processing, vol. 55, no. 7, pp. 3704–3716, July 2007.
[18] F. J. Harris, “On the use of windows for harmonic analysis with the
discrete fourier transform,” Proceedings of the IEEE, vol. 66, no. 1, pp.
51–83, Jan 1978.
[19] J. Capon, “High-resolution frequency-wavenumber spectrum analysis,”
Proceedings of the IEEE, vol. 57, no. 8, pp. 1408–1418, Aug 1969.
[20] R. Schmidt, “Multiple emitter location and signal parameter estimation,”
IEEE Trans. on Antennas and Propagation, vol. 34, no. 3, pp. 276–280,
March 1986.
[21] E. D. Di Claudio, “Asymptotically perfect wideband focusing of multiring circular arrays,” IEEE Trans. on Signal Processing, vol. 53, Part 1,
no. 10, pp. 3661–3673, Oct. 2005.
[22] E. D. Di Claudio and G. Jacovitti, “Wideband source localization by
space-time MUSIC subspace estimation,” in 2013 8th International
Symposium on Image and Signal Processing and Analysis (ISPA), Sept
2013, pp. 331–336.
[23] M. Wax and T. Kailath, “Detection of signals by Information Theoretic
Criteria,” IEEE Trans. on Acoustics Speech and Signal Processing,
vol. 33, no. 2, pp. 387–392, April 1985.
[24] A. Swindlehurst and T. Kailath, “A performance analysis of subspacebased methods in the presence of model errors: Part II - Multidimensional algorithms,” IEEE Trans. on Signal Processing, vol. 41, no. 9,
pp. 2882–2890, Sept. 1993.
[25] P. Stoica and K. C. Sharman, “Maximum likelihood methods
for direction-of-arrival estimation,” IEEE Transactions on Acoustics,
Speech, and Signal Processing, vol. 38, no. 7, pp. 1132–1143, Jul. 1990.
[26] A. Swindlehurst and T. Kailath, “A performance analysis of subspacebased methods in the presence of model errors. Part I: The MUSIC
algorithm,” IEEE Trans. on Signal Processing, vol. 40, no. 7, pp. 1758–
1774, July 1992.
[27] G. H. Niezgoda and J. Krolik, “The application of adaptive multichannel spectral estimation for broadband array processing,” in Acoustics,
Speech, and Signal Processing, 1991. ICASSP-91., 1991 International
Conference on, Apr 1991, pp. 1397–1400 vol.2.
[28] J. Li and P. Stoica, Eds., Robust adaptive beamforming, ser. Wiley Series
in Telecommunications and Signal Processing. Hoboken, NJ, USA: J.
Wiley and Sons, Oct. 2005.
[29] G. Gelli and L. Izzo, “Cyclostationarity-based coherent methods for
wideband-signal source location,” IEEE Transactions on Signal Processing, vol. 51, no. 10, pp. 2471–2482, Oct 2003.
[30] W. A. Gardner, A. Napolitano, and L. Paura, “Cyclostationarity: Half a
century of research,” Signal Processing, vol. 86, no. 4, pp. 639 – 697,
2006. [Online]. Available: http://www.sciencedirect.com/science/article/
pii/S0165168405002409
[31] E. D. Di Claudio and G. Jacovitti, “Space-time signal subspace estimation for wide-band acoustic arrays,” in 2014 22nd European Signal
Processing Conference (EUSIPCO), Sept 2014, pp. 1940–1944.
[32] O. Ledoit and M. Wolf, “A well-conditioned estimator for largedimensional covariance matrices,” Journal of Multivariate Analysis,
vol. 88, no. 2, pp. 365–411, 2004.
[33] M. L. McCloud and L. L. Scharf, “A new subspace identification
algorithm for high-resolution DOA estimation,” IEEE Transactions on
Antennas and Propagation, vol. 50, no. 10, pp. 1382–1390, Oct 2002.
[34] H. Akaike, “A new look at the statistical model identification,” Automatic
Control, IEEE Transactions on, vol. 19, no. 6, pp. 716–723, Dec 1974.
[35] K. P. Burnham and D. R. Anderson, “Multimodel inference:
Understanding AIC and BIC in model selection,” Sociological Methods
and Research, vol. 33, no. 2, pp. 261–304, 2004. [Online]. Available:
http://smr.sagepub.com/content/33/2/261.abstract
[36] G. Golub and C. V. Loan, Matrix Computations, 2nd ed. Baltimore,
USA: John Hopkins University Press, 1989.
[37] D. F. Gingras, “Robust broadband matched-field processing performance
in shallow water,” IEEE Journal of Oceanic Engineering, vol. 18, no. 3,
pp. 253–264, July 1993.
[38] R. Cicchetti, A. Faraone, E. Miozzi, R. Ravanelli, and O. Testa, “A
high-gain mushroom-shaped dielectric resonator antenna for wideband
wireless applications,” IEEE Transactions on Antennas and Propagation,
vol. 64, no. 7, pp. 2848–2861, July 2016.
[39] E. D. Di Claudio, G. Jacovitti, and A. Laurenti, “Robust direction estimation of UWB sources in Laguerre-Gauss beamspaces,” in 2011 IEEE
International Conference on Acoustics, Speech and Signal Processing
(ICASSP), May 2011, pp. 2736–2739.

15

[40] G. Su and M. Morf, “Signal subspace approach for multiple wideband emitter location,” IEEE Trans. on Acoustics Speech and Signal
Processing, vol. 31, pp. 1502–1522, Dec. 1983.
[41] M. Wax, T.-J. Shan, and T. Kailath, “Spatio-temporal spectral analysis
by eigenstructure methods,” IEEE Transactions on Acoustics, Speech,
and Signal Processing, vol. 32, no. 4, pp. 817–827, Aug 1984.
[42] S. Mohan, M. E. Lockwood, M. L. Kramer, and D. L. Jones,
“Localization of multiple acoustic sources with small arrays using
a coherence test,” The Journal of the Acoustical Society of
America, vol. 123, no. 4, pp. 2136–2147, 2008. [Online]. Available:
http://scitation.aip.org/content/asa/journal/jasa/123/4/10.1121/1.2871597
[43] T. Shan, M. Wax, and T. Kailath, “On spatial smoothing for direction-ofarrival estimation of coherent signals,” IEEE Trans.on Acoustics, Speech
and Signal Processing, vol. ASSP-33, no. 4, pp. 806–811, April 1985.
[44] S. U. Pillai and C. Burrus, Array Signal Processing. Springer Verlag,
1989.
[45] O. Pfaffel and E. Schlemm, “Limiting spectral distribution of a new
random matrix model with dependence across rows and columns,”
Linear Algebra and its Applications, vol. 436, no. 9, pp. 2966 – 2979,
2012. [Online]. Available: http://www.sciencedirect.com/science/article/
pii/S0024379511006240
[46] B. Champagne, “Adaptive eigendecomposition of data covariance matrices based on first-order perturbations,” IEEE Transactions on Signal
Processing, vol. 42, no. 10, pp. 2758–2770, Oct 1994.
[47] J. H. Manton, “Optimization algorithms exploiting unitary constraints,”
IEEE Transactions on Signal Processing, vol. 50, no. 3, pp. 635–650,
Mar 2002.
[48] X. Mestre and M. A. Lagunas, “Modified subspace algorithms for DoA
estimation with large arrays,” IEEE Transactions on Signal Processing,
vol. 56, no. 2, pp. 598–614, Feb 2008.
[49] A. Sancetta, “Sample covariance shrinkage for high dimensional dependent data,” Journal of Multivariate Analysis, vol. 99, no. 5, pp. 949 –
967, 2008.
[50] S. Kay and S. M. (Jr.), “Spectrum analysis - A modern perspective,”
Proc. IEEE, vol. 64, pp. 1380–1419, Nov. 1981.
[51] K. Hofmann-Credner and M. Stolz, “Wigner theorems for random
matrices with dependent entries: Ensembles associated to symmetric
spaces and sample covariance matrices,” Electron. Commun. Probab.,
vol. 13, pp. 401–414, 2008. [Online]. Available: http://dx.doi.org/10.
1214/ECP.v13-1395
[52] R. R. Nadakuditi and A. Edelman, “Sample eigenvalue based detection
of high-dimensional signals in white noise using relatively few samples,”
IEEE Transactions on Signal Processing, vol. 56, no. 7, pp. 2625–2638,
July 2008.
[53] C. C. Paige, “A note on a result of Sun Ji-Guang: Sensitivity
of the CS and GSV decompositions,” SIAM Journal on Numerical
Analysis, vol. 21, no. 1, pp. 186–191, 1984. [Online]. Available:
http://www.jstor.org/stable/2157055
[54] Y. Bucris, I. Cohen, and M. A. Doron, “Bayesian focusing for coherent
wideband beamforming,” IEEE Transactions on Audio, Speech, and
Language Processing, vol. 20, no. 4, pp. 1282–1296, May 2012.
[55] M. A. Doron and E. Doron, “Wavefield modeling and array processing,
part i - spatial sampling,” IEEE Trans. on Signal Processing, vol. 42,
no. 10, pp. 2549–2559, Oct. 1994.
[56] B. Friedlander and A. Weiss, “Direction finding for wide-band signals
using an interpolated array,” IEEE Trans. Signal Processing, vol. 41,
no. 4, pp. 1618–1635, April 1993.
[57] E. D. Di Claudio, “Optimal quiescent vectors for wideband ML beamforming in reverberant fields,” Signal Processing, vol. 85, no. 1, pp.
107–120, Jan. 2005.

