Super-Resolution Estimation of UWB Channels Including the Dense Component—An SBL-Inspired Approach

In this paper, we present an iterative algorithm that detects and estimates the specular components (SCs) and estimates the dense component (DC) of single-input—multiple-output (SIMO) ultra-wide-band (UWB) multipath channels. Specifically, the algorithm super-resolves the SCs in the delay-angle-of-arrival domain and estimates the parameters of a parametric model of the delay-angle power spectrum characterizing the DC. Channel noise is also estimated. In essence, the algorithm solves the problem of estimating spectral lines (the SCs) in colored noise (generated by the DC and channel noise). Its design is inspired by the sparse Bayesian learning (SBL) framework. As a result the iteration process contains a threshold condition that determines whether a candidate SC shall be retained or pruned. By relying on results from extreme-value analysis the threshold of this condition is suitably adapted to ensure a prescribed probability of detecting spurious SCs. Studies using synthetic and real channel measurement data demonstrate the virtues of the algorithm: it is able to still detect and accurately estimate SCs, even when their separation in delay and angle is down to half the Rayleigh resolution limit (RRL) of the equipment; it is robust in the sense that it tends to return no more SCs than the actual ones. Finally, the algorithm is demonstrated to outperform a state-of-the-art super-resolution channel estimator in terms of robustness in the estimation of the amplitudes of specular components closely spaced in the dispersion domain.


Super-Resolution Estimation of UWB Channels
Including the Dense Component-An SBL-Inspired Approach Stefan Grebien , Erik Leitinger , Member, IEEE, Klaus Witrisal , Member, IEEE, and Bernard H. Fleury, Senior Member, IEEE Abstract-In this paper, we present an iterative algorithm that detects and estimates the specular components (SCs) and estimates the dense component (DC) of single-input-multipleoutput (SIMO) ultra-wide-band (UWB) multipath channels.Specifically, the algorithm super-resolves the SCs in the delayangle-of-arrival domain and estimates the parameters of a parametric model of the delay-angle power spectrum characterizing the DC.Channel noise is also estimated.In essence, the algorithm solves the problem of estimating spectral lines (the SCs) in colored noise (generated by the DC and channel noise).Its design is inspired by the sparse Bayesian learning (SBL) framework.As a result the iteration process contains a threshold condition that determines whether a candidate SC shall be retained or pruned.By relying on results from extreme-value analysis the threshold of this condition is suitably adapted to ensure a prescribed probability of detecting spurious SCs.Studies using synthetic and real channel measurement data demonstrate the virtues of the algorithm: it is able to still detect and accurately estimate SCs, even when their separation in delay and angle is down to half the Rayleigh resolution limit (RRL) of the equipment; it is robust in the sense that it tends to return no more SCs than the actual ones.Finally, the algorithm is demonstrated to outperform a state-of-the-art super-resolution channel estimator in terms of robustness in the estimation of the amplitudes of specular components closely spaced in the dispersion domain.

I. INTRODUCTION
F UTURE wireless communication technologies will sup- port a variety of services with high quality requirements, Stefan Grebien is with Joanneum Research Forschungsgesellschaft mbH, 8010 Graz, Austria (e-mail: stefan.grebien@joanneum.at).Erik Leitinger and Klaus Witrisal are with the Graz University of Technology, 8010 Graz, Austria, and also with the Christian Doppler Laboratory for Location-Aware Electronic Systems, 8010 Graz, Austria (e-mail: erik.leitinger@tugraz.at;witrisal@tugraz.at).
Bernard H. Fleury is with the Institute of Telecommunications, Vienna University of Technology, 1040 Vienna, Austria (e-mail: bernard.fleury@tuwien.ac.at).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TWC.2024.3371352.
Digital Object Identifier 10.1109/TWC.2024.3371352 addressing performance metrics such as reliability, ultralow latency, high data rates, and resource-efficient use of the infrastructure [1], [2].Holistic approaches that combine different functionalities have proven to offer promising solutions to meet these requirements.Illustrative examples are integrated sensing and communications (ISAC) and radiobased simultaneous localization and mapping (SLAM) [2], [3], [4].These examples emphasize the reliance of these technologies on comprehensive, accurate channel state information.High-performance feasible parametric multi-antenna channel estimators can provide this information.

A. State of the Art
Parametric channel models typically represent multipath propagation as a linear superposition of weighted Dirac delta distributions -or spectral lines -with distinct supports in the underlying dispersion domain (delay, angle of arrival, angle of departure, Doppler frequency, and combinations thereof).Each component in the superposition is meant to represent a specular component (SC).Note that in this paper we shall use the terms SC and spectral line indiscriminately.The finite aperture of the measurement equipment imposes some limitation on the ability to resolve SCs closely spaced in the dispersion domain.
If the number of spectral lines is known, (constrained and unconstrained) maximum-likelihood (ML) methods, see e.g.[5] or subspace-based methods [6], [7] are standard super-resolution 1 tools to estimate their parameters.Expectation-maximization and related algorithms [8], [9] have proven viable approximations of the computationally prohibitive direct implementation of the constrained ML method.These estimators have in common that they do not incorporate the estimation of the number of spectral lines into the estimation problem.Schemes that perform jointly detection of the spectral lines and estimation of their parameters have been designed within a Bayesian framework [10], [11].Traditional methods combining detection and estimation select among multiple candidate models, each corresponding to a specific hypothesis on the number of spectral lines, the one that optimizes a so-called information criterion, such as the Akaike or Bayesian information criterion, and the minimum description length, see [12] and references therein.Yet, the information-based approach suffers from two shortcomings: (a) it is computationally intensive as the adopted information criterion needs to be computed first for each model candidate before a decision can be made; (b) the number of spectral lines of the selected model tends to be positively biased in non-asymptotic regimes of the signal-to-noise-ratio (SNR) and the number of observed samples [13].Hence, inference schemes designed with this approach are prone to return spurious spectral lines that have no real counterpart.Alternative penalty terms have been proposed that prevent [14] or control [15] this bias.
Model-order selection is inherently realized in sparse signal reconstruction (SSR), see [16] and references therein.SSR aims at recovering a sparse weight vector in an underdetermined linear model with a known and fixed dictionary matrix.To that end it computes an estimate of the weights as the solution to a regularized optimization problem in which the regularization term is selected to promote sparse solutions.A popular instance of SSR is basis pursuit denoising [17], also called LASSO (least absolute shrinkage and selection operator) [18], that uses an ℓ 1 -norm regularization.SSR can be formulated within the Bayesian framework as maximum-a-posteriori (MAP) estimation while imposing a sparsity promoting prior on the weight vector.Typically this prior is endowed with a hierarchical structure involving a hyperparameter for each weight.Several hierarchical models have been considered so far: gamma-Gaussian2 [20], [21], Bernoulli-Gaussian [22], [23], and generalized-gamma-powerexponential [24].This Bayesian formulation has proven to be a particularly flexible and effective tool for SSR.Since direct implementation of the estimators is typically computationally prohibitive, one has to resort to iterative schemes, often designed using variational inference methods [25], [26].
SSR can be straightforwardly applied in the context of line spectral estimation by discretizing (gridding) the dispersion domain, see e.g.[27], [28], and [29].The benefit of doing so is that the complex optimization problem needed to estimate the supports of the spectral lines is replaced by a linear programming procedure that returns a sparse estimate of the weight vector.The shortcoming is that gridding of the dispersion domain induces spectral leakage due to the resulting model mismatch.This effect can be mitigated by selecting a denser grid, yet at the cost of increasing the coherence of the dictionary matrix, which impairs the sparse reconstruction capability and increases the computational complexity.Variants of gridding methods that employ some interpolation method [30], [31], [32], [33], [34] or apply a grid refinement technique [29], [35] have been proposed to circumvent the leakage effect.
Atomic noise minimization (ATM) provides an elegant natural means to operate with a continuous, i.e. infinite, dictionary in SSR and thereby to relax the need for discretizing the dispersion domain [36], [37], [38].However, some specificities of our underlying model -namely a two-dimensional (dim.)dispersion domain and unknown colored noise -prevent a direct application of the method, see Subsection IV-A and a related discussion in [38].Moreover, numerical evidence shows that ATM requires the supports of spectral lines to be sufficiently separated in the dispersion domain in order to be able to recover them [37].In [39] an alternative is proposed that circumvents this shortcoming.
In theory, gridding-based line spectral estimation methods can be straightforwardly extended to account for continuous dispersion parameters by relaxing the discretization constraint and instead including the estimation of the support of the spectral lines in the inference process.Clearly, this approach is an instance of SSR with learning the continuous (vectorvalued) parameter of a parameterized dictionary matrix.It has been extensively pursued in connection with the Bayesian formulation of SSR [40], [41], [42], [43], [44], [45], [46].These algorithms differ in their specific design criteria, such as (i) the chosen sparsity-inducing hierarchical prior model, e.g.gamma-Gaussian [40], [41], [42], [43], [44], Bernoulli-Gaussian [45], [46], (ii) the assumed absence [40], [46] or presence [41], [44] of correlation among the weights of the spectral lines, and (iii) whether point estimates [40], [42], [43], [44], [45] or posterior probability density functions (PDFs) of the dispersion parameters of the SCs are inferred [46].Experimental evidence shows that the algorithms computing point estimates of the supports of spectral lines show a positive bias in the number of detected spectral lines, i.e. are prone to detect spurious spectral lines.Including inference of the posterior PDF of the supports allows for mitigating this bias, yet at the cost of an increased computational complexity [46].We remark that (iterative) SSR methods that apply grid refinement techniques such as [29], [35] can be viewed as particular instances of SSR methods with continuous-parameter learning, which adapt their inherent restricted range of the dictionary parameter during the iterations. 3 The above SSR methods with continuous-parameter dictionary learning include an inherent pruning procedure that determines which ones among the columns of the dictionary matrix are inferred as relevant and switch the others off, see e.g.[47] and [48].It is shown in [44] that the number of detected spurious spectral lines can be significantly reduced by suitably adapting the threshold of the pruning stage.The analysis provided there relies on some heuristic, yet realistic, assumptions that allow for approximating the probability of detecting a spurious line with the probability that the maximum of a continuous χ 2 random field exceeds the selected threshold [49], [50].The analysis shows that a prescribed probability of detecting spurious lines can be guaranteed, provided the threshold increases as C + log n + 1 2 log log n where n is the number of observation samples and C is a constant that depends on that probability [14], [44].Numerical analyses have shown that using this adapted threshold leads to almost vanishing bias in the number of detected SCs in medium and high SNR regimes with a tendency to underestimate said number in the low SNR regime, see also Section VII.
In recent years, an extension of the channel model has been considered, that includes a DC [9], [51], [52].The DC incorporates diffuse scattering as well as SCs that cannot be resolved with the finite aperture of the measurement equipment.Including the estimation of the DC can improve the accuracy of the estimation of the parameters of resolved SCs [9].

B. Contributions of the Paper
We propose an iterative algorithm that performs combined detection and estimation of SCs and estimation of the DC plus additive white Gaussian noise (WGN) (AWGN) in singleinput-multiple-output (SIMO) ultra-wide-band (UWB) multipath channels. 4The algorithm resolves the SCs in the delay-angle-of-arrival (angle) domain.The contributions of this paper are as follows: • We model the impact of the DC and AWGN as colored noise, so that the problem becomes that of line spectral estimation [54] in such noise when the relative delays that the (UWB) complex envelope of the sounding wave exhibits when it is sensed by the elements of the antenna array cannot be neglected.• The design of the algorithm is inspired by the SBL approach [20].The probabilistic model is extended by assuming that the weights of the spectral lines are independent circularly-symmetric complex Gaussian random variables with unknown variances.In a first stage ML estimation of the variances and all other parameters but the weights is performed after integrating out said weights.These estimates are then used to compute a tractable (Gaussian) approximation of the weights' posterior PDF.The algorithm computes these two stage, the former one in an iterative fashion.• We suitably modify the threshold inherent to the above ML estimation stage to meet a prescribed probability of detecting spurious spectral lines.To do so we apply results from extreme value analysis [49], [50].• Using synthetically generated observation data we study in-depth the behavior of the proposed algorithm and especially how the adapted threshold affects its performance.• We compare the performance of the proposed algorithm with a maximum-likelihood estimation (MLE) algorithm inspired by [9] but with different scheduling and adapted threshold according to [14].• We apply the algorithm to UWB measurement data collected in an indoor environment.A simple ray-tracing tool is used to identify plausible propagation paths that can be associated to the SCs detected by the algorithm.The remainder of the paper is organized as follows: In Section II we present the generative signal model for the considered SIMO measurement set-up.Section III describes the probabilistic signal model for inference.We derive the proposed algorithm in Section IV.Section V addresses the analytical correspondence between the probability of detecting spurious spectral lines and the threshold of the ML estimation stage.Section VII reports results from numerical and experimental studies.Concluding remarks are provided in Section VIII.

A. Continous-Time Signal Model
The experimental measurement setup consists of an UWB transceiver operating in an indoor environment.The transmitter (Tx) is equipped with a single antenna, while an antenna array with colocated elements is emulated at the receiver (Rx) using a single antenna mounted on a positioning table.
Signals are represented by means of their complex envelope with respect to a center frequency f c .For the sake of simplicity we assume horizontal-only propagation. 5Under the planewave assumption, provided that the location p of the Rx antenna is restricted to a region R sufficiently confined around a reference point p ∈ R, the antenna output signal can be expressed as In this expression t is the time variable, while where s(t) is the transmitted signal with bandwidth B and g φ, p = [cos(φ) sin(φ)](p − p)/c with c denoting the speed of light, expresses for a plane wave incident with angle φ ∈ [−π, +π) the wave's excess (propagation) delay at p relative to its delay at the reference point p.The function h(τ, φ) ∈ C defined on R × [−π, +π) characterizes the spread in (relative) delay τ and angle φ of the signal sensed at p. Finally, w(t; p) is a spatially and temporally white Gaussian noise field with double-sided power spectral density N 0 /2.Clearly, both h(τ, φ) and z t; τ, φ, p and thereby r(t; p) given by (1) also depend on the reference point p.Since this point is kept fixed, the used notation discards this dependency for the sake of shortness.Sufficient conditions for the right-hand expression in (1) to accurately represent r(t; p) are that (a) the region R is far away enough from the Tx and the objects in the environment with which the dominant multipath components interact along their path, such as walls, boards, etc., and (b) that the spread function h(τ, φ) stays constant over the bandwidth (frequency aperture) of the sounding signal.The latter assumption implies that the electromagnetic properties of said objects, like reflection and transmission coefficients, are nearly constant over the sounding bandwidth.Note that the directivity of the Rx antenna is absorbed in h(τ, φ).
The rationale behind the selection of model ( 3) is as follows.The SCs originate from electromagnetic interactions with objects in the environment that are essentially nondispersive, such as line-of-sight (LOS) propagation, specular reflection and transmission, and can be resolved with the used aperture.The DC incorporates the contributions from all other interactions, e.g.diffuse scattering and diffraction.It also includes components from specular interactions that cannot be resolved with the used aperture.
The antenna array is emulated by measuring r(t; p) in (1) at M positions p m ∈ R, m ∈ {1, . . ., M } ≜ M, see Fig. 1.These positions specify the array layout.Without loss of generality the center of gravity of the layout is set to coincide with the reference location, i.e. p = M −1 M m=1 p m .The angle o depicted in the figure determines its orientation.
Inserting the decomposition (3) in (1) yields for the antenna output signal at the mth location We have defined r m (t) = r(t; p m ) and w m (t) = w(t; p m ) for short.

B. Discrete-Frequency Signal Model
The signals r m (t), m ∈ M are Nyquist filtered, Fourier transformed, and then synchronously and uniformly sampled with frequency spacing ∆ over the bandwidth B to collect for each branch m N = B/∆ samples that are arranged in a N -dim.vector y m .These M vectors are then stacked to form the N M -vector y (5) With Z(f ; τ, φ, p m ) denoting the Fourier transform of z(t; τ, φ, p m ), i.e., where S(f ) is the Fourier spectrum of s(t), the mth entry in the vector in (6) reads i.e., it contains the Fourier-transformed samples collected at the mth antenna element.The N M -vector n = v + w in ( 5) aggregates the vectors v and w that collect the samples (arranged in the right order) corresponding to, respectively, the integral term and the noise term in (4) when m ranges in M. From the assumptions on the DC, v is a complex circular symmetric M N -dim.Gaussian random vector with zero mean and with (m, m ′ ) ∈ M 2 .From the assumptions on the noise measurement process w is a complex circular symmetric Gaussian random vector with covariance matrix σ 2 I N M where σ 2 = N 0 /T s and I (•) is the identity matrix of dimensions specified by the number given in the subscript.We assume that v and w are uncorrelated.As a result n is a circularly symmetric Gaussian random vector with covariance matrix

C. Selected Model for the DC
We impose some structure on the covariance matrix Q in (10) via some assumptions on the behaviour of the DC and simplifying approximations in the derivations of the submatrices in (9).This structure will ensure the feasibility of the estimation algorithm.
b) In the computation of ( 9) we discard the second occurrence of the term g φ, p m in (7), i.e., Z(f ; τ, φ, p m ) = Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.e j2πfcg(φ,p m ) S(f )e −j2πf τ .This step amounts to adopting a narrowband representation that neglects the relative delays across array elements of the modulating signals of incident waves [9], [58].It follows from Assumptions a) and b) that the covariance matrix of v factorizes as with ⊗ denoting the Kronecker product [9], [58]. 6The first factor is the spatial correlation matrix being the steering vector of the array, and more generally the spatial steering vector of the sounding equipment.The second factor is the delay correlation matrix with a f (τ ) = S i∆ e j2πi∆τ : i = −(N − 1)/2, . . ., (N − 1)/2 T ∈ C N ×1 denoting the frequency steering vector of the sounding equipment.c) Experimental evidence shows that the DPS typically exhibits an exponentially decaying tail [9], [55], [59] and a smooth onset [55], [60].This behavior is well represented by a truncated and normalized gamma PDF given by where u(τ ) is the unit step function, Γ(•) is the gamma function, and ϑ = [β θ ξ] collects the onset (β > 0), scale (θ > 0), and shape (ξ > 0) parameters.The normalization constant a > 0 guarantees that p(τ ; ϑ)dτ = 1.The range of the parameter β is restricted in such a way to ensure that the integral of the truncated tail of the gamma PDF is negligibly small, i.e., a(β) ≈ 1 for any such values of β. d) We neglect the spatial correlation across antenna elements, i.e., we set Q s = I M [9], [61].This choice provides a good approximation of Q s under the assumption of uniform APS for the antenna-element spacings used in practice. 7y combining Assumptions a)-d), the covariance matrix R takes the following form: where η = [σ 2 P ϑ] with ϑ defined above and Q f given in (13) with p(τ ) = p(τ ; ϑ) according to (14).Hence, because

III. SPARSE BAYESIAN FORMULATION
If the number of components K of the model ( 5) were known, the vectors of dispersion parameters ψ and complex amplitudes α of the SCs and the parameter vector η of colored noise could be inferred using a standard MAP or ML estimation technique.Since we can view the family {Z( ψ)} ψ∈Ψ as a continuous dictionary, atomic-norm methods seem at first glance to be an inference method particularly tailored to our model.However, as detailed in the discussion of Section IV, some specificities of the model prevent a direct application of these methods.
We propose an approach inspired from SBL [20] for SSR to include the estimation of the unknown K.This approach requires a two-fold modification of the generative signal model that we address below.

A. Discrete-Frequency Signal Model for Inference
In a first step the initial generative signal model ( 5) is modified as follows.The number of hypothetical SCs is set to a fixed number, say L. Parameter L is a design parameter that is selected large enough so that K ≤ L. In addition, L ≪ N M .Actually we only need that L ≤ M N .The further restriction ≪ is for feasibility issues.Similarly as in Subsection II-B, we define the vector With these modifications, we arrive at the discrete-frequency signal model given by where with z(ψ l ) defined similarly to (6).Under the assumptions made in Subsection II, the likelihood function of this model reads with det(•) denoting the determinant of a matrix.The second step consists in specifying a hierarchical prior for each entry α l in form of a Gaussian scale mixture.Specifically, we define where Note that all entries in γ have the same prior with PDF f (γ).These entries and their prior are referred to as hyperparameters and hyperprior, respectively.We postulate priors for the parameter vectors ψ and η with respective PDFs f (ψ) and f (η).With these specifications, the probabilistic model for inference reads Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Inference Method
The proposed method is inspired by SBL [20].First it computes a MAP estimate of ψ, η, and γ from the joint posterior of these random vectors; then it uses these estimates to infer an approximation of the posterior distribution of α.
The posterior PDF f (ψ, η, γ|y) is obtained from (20) by marginalizing out the complex amplitude vector α, i.e., where with and The MAP estimates of ψ, η and γ are then computed using (21).From (20) we get which is readily shown to be Gaussian with mean and covariance matrix The approximate posterior PDF of α results by plugging the MAP estimates of ψ, η and γ in (24), and thus in (25) and (26).
In our design, we select non-informative improper priors for ψ, γ and η: With this selection, the above MAP estimates coincide with the ML estimates and the posterior PDF of α is inferred using the approximation f (α|y, ψML , ηML , γML ).

IV. ITERATIVE DESIGN OF THE ESTIMATOR
Since the ML estimator in ( 27) cannot be calculated analytically, even though the likelihood function is given in an analytical form, and a direct numerical solution is computationally prohibitive, we resort to a sequential update of the parameter vectors ψ, η, and γ resulting in the estimates ψ, η, and γ.
1) Estimation of the Supports of the Spectral Lines: Inserting the current estimates η, and γ in (27) the new estimate of ψ is computed to be 2) Estimation of the Parameters of Colored Noise: Similarly, the new estimate of η is computed based on the current estimates ψ, and γ to be 3) Estimation of the Hyperparameters: Finally, given the current estimates ψ and η, the new estimate of γ is updated according to In the sequel we consider instead of (30) a sequential method in which the estimate of each entry in γ is updated while the estimate of the other entries are kept fixed [47]: with κ = 1.In this expression with where diag(•) describes a square diagonal matrix with the elements of the vector given as an argument on the main diagonal and ψl = [ ψ1 Note that the computation step of γl (31) contains a condition that determines when the lth spectral line shall be discarded (γ l = ∞).
4) Estimation of the Weights: Inserting the estimates ψ, η, and γ in (24) yields the Gaussian PDF with mean (see (25)) and covariance matrix (see (26)) that is used as an approximation of the posterior PDF of α.
In (39), 5) Fitting of the Pruning Threshold κ: Numerical experiments have shown that the iterative algorithm obtained in the above subsections overestimates the number of spectral lines and thereby returns estimates of spurious components.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
This bias in the number of detected components increases when either the SNR or the number of samples increases [40, Subsec.V.A], [44].Following the approach adopted in [14] and [62], we increase the initial threshold κ = 1 in the pruning condition (31) to κ = κ * > 1.The value κ * is set in such a way to reduce the bias.The next section describes in detail this procedure, which yields κ * given in (48).

A. Discussion
The updating step (28) in its form looks very similar to the classical unconstrained (also called stochastic) ML estimator in sensor array signal processing [5] with the additional assumption that the precision matrix of the weights be diagonal, i.e. equal to Γ as a result of the gamma-Gaussian hierarchical model. 8Despite the resemblance (28) is not an instance of unconstrained ML estimation.Unconstrained ML estimation requires a scenario where at least as many observations (snapshots, assumed uncorrelated) as the number of sensors are collected, 9 while in our scenario only one observation, i.e. y in (17), is available.Our estimator also deviates from being an instance of SBL [20] in three respects: (a) the underlying model of SBL is undetermined, which is not the case for our model (17) with L ≪ N M ; (b) the "dictionary matrix", namely S(ψ) in (17), is not fixed but is parameterized by the continuous parameter vector ψ that is estimated; and (c) the inherent threshold of SBL is adapted to control the probability of detecting spurious SCs. 10 Our method belongs to the class of parametric estimators in the nomenclature introduced in [27].
Atomic noise minimization (ATM) provides an elegant, natural means to operate with a continuous, i.e. infinite, dictionary in SSR [36], [37].At first glance this method looks promising for dealing with the continuous dictionary {S(ψ)} ψ∈Ψ in our problem at hand.However, some specificities of the generic model (5) prevent its straightforward application to our scenario.Note that ATM primarily "denoises" the observed signal with the estimation of the spectral lines being subsequently performed based on this denoised signal.While the estimation problem can be solved with an exact semi-definite program when the dispersion domain is one-dim.[36], only an approximate such program could be formulated to date for higher dim.dispersion domains [38]. 11In addition, ATM operates on Nyquist-sampled signals and requires knowledge of the noise characteristics, e.g. its spectral height when noise 8 Sensor array signal processing considers a signal model similar to (17) where the entries of y are the outputs of an array of sensors, Z(ψ) is the array response matrix, ψ and α contain respectively the dispersion parameters and the amplitudes of the sources, and n is the measurement noise vector [5].The number of sources is assumed to be known and smaller than the number of sensors in order for the model parameters to be identifiable.In practice the number of sources is estimated using an additional model-order selection procedure based on an information theoretic criterion, see Section I. 9 This condition ensures that the sample array covariance matrix has full rank, which is a mandatory condition in the derivation of the unconstrained ML estimator. 10Strictly speaking, SBL is derived under the assumption of AWGN.It can be straightforwardly applied when noise is non-white, by merely applying a whitening filter first. 11In [38] the matrix-enhancement-matrix-pencil method [63] is used to compute estimates of the support of spectral lines from the denoised signal.is white.These conditions do not hold in our application scenario.

V. COMPUTATION OF THE PRUNING THRESHOLD
To compute the threshold value κ * we adapt the approach described in [44] to our application scenario; see also [14] for a similar approach applied to constrained ML estimation.To make it tractable the analysis is carried out under the following assumptions.
Assumption 1: The spatial and frequency apertures [64] of the sounding equipment are centro-symmetric12 [65].Furthermore, a f = a f (0), see text below (13), fulfils J a f = a * f , where J is the exchange or reversal matrix [54,Sec. 4.8].The covariance matrix R in ( 15) is known.
It is shown in [53] that as a result of the first part in the assumption the matrix Q in ( 11) is centro-hermitian 13 [65] and therefore R in (15) too.The next assumption reflects an empirical evidence based on extensive simulations of the proposed algorithm.
Assumption 2: Asymptotically as the dimension M N grows large the estimator in Section IV with κ = 1 exhibits the following behaviour: (a) it resolves all K active SCs and accurately estimates their parameters, i.e. without loss of generality ψl ≈ ψl for l = 1, . . ., K; (b) it computes estimates ψl , l = K + 1, . . ., L of L − K (spurious) SC components in such a way that with high probability z( ψl ) is nearly orthogonal to any columns of Z( ψl) for each l = K + 1, . . ., L.
As a result of Assumption 2, as M N grows large, with high probability (32) and ( 33) can be approximated for l = K + 1, . . ., L as ζ l ≈ ζ( ψl ) and ρ l ≈ ρ( ψl ), respectively, where we have defined Therefore, the probability that the algorithm decides that the lth component (l = K + 1, . . ., L) is active, i.e. γl < ∞ in (31), with threshold set to κ is close to when M N is sufficiently large.Let us consider the circularly-symmetric complex Gaussian random field on Ψ defined as with ψ • ∈ Ψ.Then, (40) can be recast as Theorem 1: Under Assumption 1 we have the asymptotic equivalence Furthermore, Here, Λ(ψ • ) is the non-negative definite matrix given in (56) and with R given in (16) and ȧf (τ ) = ∂a f (τ )/∂τ .The term incorporates the impact of the array aperture, while b 1 (τ ) and b 2 (τ ) account for the impact of the frequency aperture (the spectrum S(f )) and colored noise.
Proof: See Appendix A.
The function of κ in the asymptotic equivalence (43) provides a tight approximation of P sp (κ) versus κ for κ sufficiently large.Thus, by taking the inverse of that function and evaluating it at a target probability value, say ϵ, we obtain a threshold, say κ ⋆ = κ ⋆ (ϵ), that yields P sp (κ ⋆ ) close to ϵ.The next lemma essentially gives this inverse function.

A. Examples
We illustrate the right-hand expression in (44) with two examples.We consider a scenario with a sounding signal exhibiting a constant spectrum over its bandwidth, i.e., Assumption 1 is fulfilled, and AWGN only.In this case (44) becomes Proof: In this case, b 1 (τ ) = 1 and b(τ □ If furthermore the array is uniform, square, of dimensions M ′ × M ′ , and with equal inter-element spacing w > 0, the above expression further simplifies to Proof: The square aperture function of a rectangular uniform array is given as □ In [67] classical detection theory for detecting a single spectral line with unknown frequency embedded in AWGN is applied to obtain a threshold that is derived by analyzing the maximum of the samples of a periodogram taken on the Fourier lattice [68].In [53] it was shown that such a classical threshold results in a positive bias in the number of detected spurious components that is avoided when the threshold given in (48) is employed.In [69], an improved threshold is proposed that reduces the positive bias in the number of spurious components too.

VI. IMPLEMENTATION A. Algorithmic Routine
The pseudocode of the proposed algorithm is given in Algorithm 1.It has two main stages: a search and a refine procedure, described in Procedure 1 and Procedure 2, respectively.After initialization the two procedures are executed sequentially in a do-while loop until a stopping criterion is met.Specifically, Algorithm 1 implements a bottom-up strategy: starting with an empty model, i.e.L = 0, at each iteration of the do-while-loop Procedure 1 searches and adds a candidate SC in the current pool of so far buffered candidate SCs.Procedure 2 estimates and/or re-estimates the parameters of all candidate SCs in the pool, and possibly removes candidate SCs to finally yield an updated pool of L candidate SCs.The algorithm terminates once the number L of SCs in the pool and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 2. Comparison of the probalilities Psp(κ) and P mi (κ) with their respective relative frequencies computed from 1000 trials with SNR as a parameter.The three scenarios corresponding each to specific assumptions on Q are described in the text.

Algorithm 1 Main
Input : Measurement vector y Output: L, ψ, η, μ, and  32) and ( 33) append ψl to ψ and γl to γ The initial iterations in the do-while-loop of Algorithm 1 are executed while considering measurement noise only, i.e. by using R in (15) with P set to zero whenever R occurs in the update equations of Procedures 1 and 2. This is carried out until a first candidate SC in the pool is pruned in Procedure 2 or L reaches a predefined maximum number L, in which case Procedure 3 is executed to initialize the parameters of the DC.Once the initialization is completed, the noise variance estimate σ2 aggregates a contribution from the DC.The total estimated power over the bandwidth computed from this value is distributed evenly between noise and the DC in Procedure 3.This explains the factor 1/2 in Lines 3 and 5. From then on

B. Complexity Analysis
The operations that mainly affect the computational complexity of Algorithm 1 are, respectively, the updating step for the dispersion parameter vector (28) and the updating step for the noise parameter vector (29) in Procedure 2. By making use of the matrix inversion and matrix determinant lemmas and exploiting the Kronecker structure of the covariance matrix in (15), the computational complexity of the former updating step is O((N M ) 2 L), while that of the latter is O((N M ) 2 L + N 3 ) (whether the first or the second term is dominating the other depends on the number of frequency samples N versus the number of antennas M times the number of estimated components L).Note that the bottom-up approach steadily increases L, starting from an empty model.Thus, the approach avoids unnecessary computation.The computational complexity of the modified MLE algorithm used as a benchmark in Section VII-A4 is O((N ) 2 M L) for updating the dispersion parameter vector and O(N 3 ) for updating the noise parameter vector.Note that the RiMAX algorithm in [9] shares the same complexity.

VII. NUMERICAL AND EXPERIMENTAL RESULTS
To validate the proposed algorithm, we first test it in Subsection VII-A with synthetically generated measurements according to the model in ( 5) with a covariance matrix R given in ( 9) and (15).Then, in Subsection VII-B we apply the algorithm to measurements collected in an indoor environment.

A. Synthetic Radio Channels
In this study, the signal spectrum S(f ) has a root-raised cosine shape with roll-off factor 0.6 and bandwidth B = 1.6 GHz centered at f c = 6 GHz.The 3-dB bandwidth B is 1 GHz and yields the RRL 1/ B = 1 ns.Each numerical investigation involves 1000 simulation trials.In each trial the Gaussian DC vector v (see text below ( 8)) is generated using (14) with β = 1m/c, θ = 5 ns and ξ = 1.8.The power P is specified through the SC-to-DC-ratio SDR = 10 log 10 1 M ∥ k∈K αk z( ψk )∥ 2 /(P B) .In addition, the Gaussian noise vector w is generated with component variance σ 2 specified through the signal-to-noise ratio SNR = 10 log 10 ( 1 M ∥ k∈K αk z( ψk )∥ 2 + P B)/σ 2 .Parameter settings common to all experiments considered in the paper are reported in Table I.
1) Empirical Substantiation of the "Near-orthogonality" Assumption: This study presents empirical evidence supporting Assumption 2. We consider a synthetic channel with a single SC, i.e., (5) with K = 1.The settings specific to the experiment are summarized in Table II.The dispersion parameters of the SC are selected as follows: its delay is fixed to τ1 = τ = 10 ns and its angle φ1 = φ is drawn uniformly over [0, 2π) for each trial and independently across trials.The respective powers of noise, the DC and the SC are set such that SDR = −5 dB and SNR = 20 dB and R is computed using (15).
To substantiate the "near-orthogonality" property claimed in Assumption 2 in each trial the cross-correlation coefficient is calculated for any pair (i, j) of indices of SCs detected by the algorithm and the mean and variance of these figures are obtained.The latter quantities are averaged over the 1000 trials to yield 0.0123 and 0.0401, respectively.As a note the average number of detected SCs is 9.52.We also plotted (not reported here due to space constraints) the estimated dispersion vectors of the detected SCs in their domain Ψ.By visual inspection we could qualitatively observe that vectors located outside an elliptically shaped region centered at the dispersion vector [τ φ] of the active SC look uniformly distributed.The main axes of the boundary ellipse are set equal to 5 times the root-Cramér-Rao lower bounds (CRLBs) for the estimation of the delay and angle.This choice ensures that estimated dispersion vectors located outside the elliptical region are very unlikely (specifically with probability equal to 0.0001) to be a noisy estimate of [τ φ].
2) Detection and Estimation of a Single SC: In this study, we first validate empirically the expression in (43) as an approximation of the probability of detecting spurious SCs as well as an expression approximating the probability of not detecting an active SC that we introduce now.In a single-SC scenario, the probability of missed detection can be approximated in the asymptotic regime N M → ∞ by where ψ is the estimated dispersion vector of the detected SC.In this regime the distribution of 2|x( ψ)| 2 can be approximated by a non-central χ 2 distribution with 2 degrees of freedom and non-centrality [44].Making use of this result, the probability of missed detection in a single-SC scenario is approximated by [44] To numerically assess the accuracy of using ( 42) and ( 54) as approximations of the probabilities of detecting spurious SCs and missing an active SC, respectively, we modify the settings as summarized in Table II.The threshold κ * of our algorithm is a varying parameter.Other not explicitly mentioned settings stay as described in Subsection VII-A1.Note that by keeping the delay of the SC fixed the non-centrality parameter stays constant and equal to η = {8.2,12.2, 16.7} dB corresponding to the three SNR values.
Fig. 2 shows a comparison of P sp (κ) in (42) (dashed lines) and P mi (κ) in (54) (dash dotted lines) with the relative frequencies of, respectively, detecting a spurious SC (solid lines) and missing the active SC (dotted lines) computed from 1000 trials.To compute the latter quantities we count the occurrence of two events that we now define.First we specify a rectangular region in Ψ centered at the dispersion vector of the active SC and with sides equal to 5 times the square root of the respective CRLBs [70].The event "false detection" occurs if the estimated dispersion vector of at least one detected SC lies outside the region.The event "missed detection" occurs if no SC is detected or the estimated dispersion parameters of all detected SCs lie outside the region.The study is conducted under two assumptions on the covariance matrix R in (10) used in the generative model: R has the simplified form (15) (Fig. 2a and Fig. 2b) and R has the general Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II
EXPERIMENT SPECIFIC SYSTEM PARAMETERS Fig. 3. Wideband versus narrowband detection and estimation of SCs in the AWGN channel: (a) mean number of detected SCs, (b) root mean-square error (RMSE) of the distance estimates of detected SCs, and (c) RMSE of their angle estimates obtained with the proposed algorithm (in blue with crosses) and its narrowband version (in red with pluses) as a function of the array size M .form (10) (Fig. 2c). 14Furthermore, under the first assumption we distinguish between the two cases where the matrix R is known (Fig. 2a) and unknown (Fig. 2b) to the algorithm, and thus is estimated in the latter case.
We see in Fig. 2a and Fig. 2b that when the matrix R selected in the generative model matches that used in the design of the algorithm, whether the algorithm knows or does not know said matrix has little impact on its performance.Fig. 2c shows that when there is a mismatch, it only marginally affects the performance of the algorithm.Specifically, a comparison with Fig. 2b shows that at large SNR values the number of spurious SCs is slightly increased due to the mismatch.
3) Wideband Versus Narrowband Detection and Estimation of SCs in the AWGN Channel: State-of-the-art detection and estimation schemes are traditionally designed based on the narrow-band assumption, which neglects the second occurrence of g(φ, p m ) in ( 7) [8], [9].In this study we show that neglecting this term in the proposed algorithm leads to an increase of the number of detected spurious SCs as the size of the array increases.To quantitatively assess this effect, we modify the simulation scenario as summarized in Table II; σ 2 is assumed known.
In this study and the subsequent ones we adopt the widely used convention in the radar community to convert (propagation) delays in their corresponding equivalent (propagation) distances.Fig. 3 depicts results obtained from 1000 simulation trials that illustrate the behavior of the proposed algorithm (in blue with crosses) and of a simplified (narrowband) version of it that neglects the second occurrence of g(φ, p m ) in (7) (in red with pluses) as a function of the array size M .Fig. 3a, Fig. 3b and 3c report respectively the mean number of detected SCs, the RMSE of the distance estimates, and the RMSE of the angle estimates.When M is increased beyond 5, the narrowband assumption is violated and the narrowband version of the algorithm detects additional spurious SCs with dispersion vectors located in the vicinity of that of the active SC.By contrast, the RMSEs achieved with the proposed algorithm decrease slightly as M is increased. 15) High-Resolution Capability of the Proposed Algorithm: To study the super-resolution capability of the proposed algorithm, we consider a scenario with K = 2 SCs with a controlled separation of their respective dispersion vectors.The parameter vector ψ1 of the first SC is drawn uniformly over Ψ.The parameter vector of the second SC is then set to either ψ2 = ψ1 + [∆d/c 0] or ψ2 = ψ1 + [0 ∆φ].The spacings ∆d and ∆φ are fractions of the RRL in, respectively, distance (c/ B ≈ 0.3 m) and angle (56 • ) [71].The complex amplitudes of both SCs have unit magnitude and their respective phases are drawn uniformly and independently.Other system parameters are summarized in Tables I and II.
In Fig. 4, we compare the performance of the proposed algorithm (circles) with that of a MLE algorithm (trian-  gles) inspired by [9].The MLE algorithm that we use as a benchmark differs from that in [9] in two respects: the scheduling and the thresholding that both coincide with those implemented in the proposed algorithm.As mentioned in the manuscript, the method used to obtain the threshold is inspired from that in [14]. 16Performance versus spacings ∆d and ∆φ are depicted as, respectively, blue solid and red dashed curves.
The first two columns of Fig. 4   capability of the algorithms is restricted by the SDR.Further worth mentioning is that both algorithms tend to underestimate the number (K = 2) of active SCs when the spacing is reduced.
Columns three and four of Fig. 4 depict the RMSE of, respectively, the distance and angle estimates, provided exactly two SCs are detected.We associate the two detected SCs with the true SCs by means of the optimal subpattern assignment (OSPA) metric [72].To be able to use the metric we normalize the estimated distances and angles with the RRL in distance and in angle, respectively.The root of the sum of the CRLBs (delay and angle) of the two SCs are also depicted (lines with stars).The estimates returned by both algorithms approach their respective root CRLB, provided the spacing of the two SCs in Ψ is large enough.
The last column in Fig. 4 presents the mean of the absolute value of the complex amplitudes ⟨ l |α l |⟩ (empty markers) and the RMSE of the absolute value of the complex amplitudes (filled markers) again provided exactly two SCs are detected.The two algorithms perform similarly in estimating the dispersion parameters; however our algorithm outperforms the modified MLE algorithm in estimating the complex amplitudes when the two SCs are closely spaced.This distinct behavior results from the specific structures of the algorithms: the modified MLE algorithm computes a least-squares estimate of the amplitudes, while the proposed algorithm computes a linear minimum mean-square error (MMSE) estimate based on the hyperparameter estimates, see (38) and (39).When the estimated dispersion vectors of the two detected SCs are closely spaced, the least-squares estimator computes the inverse of an ill-conditioned matrix, while the linear MMSE estimator regularizes this matrix.

B. Measured Radio Channels
For the experimental study we used a channel sounding equipment that transmits an m-sequence of 7 GHz bandwidth at 6.95 GHz carrier frequency.Details about the equipment can be found in [73].After applying standard pre-processing steps (subtracting the cross-talk and equalizing the system response), the resulting signal is input to a filter with a root-raised-  II for the parameter settings.
The room where the measurements were performed is depicted in Fig. 5. Based on a layout of it, see Fig. 7 in Appenix D, a classical mirror source method [74] computes the positions of predicted virtual sources associated with rays from the Tx antenna to the center of gravity of the Rx array positions that undergo up to 5 reflections on walls or large objects (windows, boards).To each such predicted source corresponds a predicted SC with dispersion vector computed from the position of the source, see Appendix D.
The following analysis concerns measurements obtained with Rx position p 1 depicted in Fig. 7 in Appendix D. Fig. 6(a) depicts the estimated DAPS computed from the received signal [74].Note that this power spectrum incorporates the smoothing function of the aperture of the measurement equipment [64].This will be the case for all power spectra considered in this study.The red crosses and blue diamonds mark the estimated dispersion vectors of the SCs detected by the algorithm and the predicted SCs, respectively.To each detected SC we associate (possibly no, one, or more than one) predicted SC as follows.A predicted SC is associated to a detected SC if their respective distances and angles are no more apart than, respectively, 10 cm (1/3 of the RRL) and 5 • (1/10 of the RRL), see Appendix D for the rationale behind this choice.To most of the detected SCs, a unique predicted SC is associated in this way.The two detected SCs with dispersion parameters 35 m and 80 • are associated with the same predicted SC.An association could not be made for four detected SCs. Fig. 6(c) shows the estimated DAPS computed from the residual signal y−Z( Ψ) μ.It also includes the dispersion vectors of detected SCs (blue diamonds) that could not be associated with any detected SC. Clearly, the strong peaks in the estimated DAPS depicted in Fig. 6(a) have vanished.Fig. 6(b) shows the estimated DPS computed from the original signal (solid blue with crosses) and from the residual signal (solid red with pluses), as well as the theoretical DPS of the DC given in ( 14) with η = η plus the estimated noise variance σ (solid black with stars).The first two DPS result from averaging the respective DAPS over the angle domain.The DPS obtained from the residual and reconstructed signals match well.This empirically justifies our choice of model (14).Finally, Fig. 6(d) depicts the estimated APS computed in a similar way as the DPS depicted in Fig. 6(b).The first two estimated APS are obtained by averaging the respective DAPS over the delay domain.It can be seen that the estimated APS of the residual signal is nearly constant over the angle domain.

VIII. CONCLUSION
In this paper, we derive and analyze a super-resolution algorithm for detecting and estimating specular components as well as estimating the power spectrum of the dense component plus noise in an ultra-wide band single-input-multiple-output (SIMO) multipath channel.Estimated parameters are among others the delay, angle-of-arrival, and complex amplitude of the detected specular components as well as the parameters of a parametric model of the delay power spectrum characterizing the dense component.The design of the algorithm is inspired by sparse Bayesian learning.As a result it embodies a pruning condition that determines whether a candidate specular component is considered active or not.The threshold of the pruning condition is adapted to control the probability of detecting spurious specular components.
Numerical studies in a synthetic environment show that the simplifying assumptions underlying the derivation of the algorithm are realistic and that the relative frequencies of detecting spurious specular components and missing active specular components are close to the respective probabilities derived theoretically.These studies also demonstrate several virtues of the algorithm: (a) its ability to still detect and accurately estimate specular components, even when their separation in delay and azimuth is down to half the Rayleigh resolution limit of the equipment; (b) it is robust in the sense that it tends to detect no more specular components than the actual ones.An experimental study illustrates the ability of the proposed algorithm to accurately infer the dispersive characteristics (in delay and angle of arrival) of the UWB SIMO channel.Owing to this high efficiency the proposed algorithm has promising potential applications in all areas of wireless communications that exploit extended channel state information, such as integrated sensing and communications (ISAC) and radio-based localization.

APPENDIX A PROOF OF THEOREM 1
Proof: As shown in [53] it follows from Assumption 1 that the real and imaginary parts of the Gaussian field x(ψ • ) in ( 41) exhibits the following properties: 1) They have equal constant variance: 2) They are independent: is a random χ 2 field on Ψ with two degrees of freedom [49], [50]. 17Note that unless the DC vanishes, i.e.P = 0, see (15), the Gaussian field x(ψ • ) is non-stationary and so is 2|x(ψ • )| 2 .The probability that the χ 2 field exceeds a threshold is asymptotically equivalent to the probability of the field's excursion above the threshold when said threshold grows large [49], [50].Specifically, by applying Weyl's tube formula [50, Theorem 3.
In this expression As shown in Appendix C, the entries in ( 56) read The right-hand side in (44) follows then from ( 57)-( 59).□
1) Second-Order Partial Derivatives of x(ψ) w.r.t.τ : One can easily check that because of (61) the second and third terms of (60) vanish in this case.As a result, we can write .
The above expression is a delay-dependent loss factor that depends of the structure of the noise vector n.Note that if n is white, b 2 (τ ) = 1, otherwise b 2 (τ ) < 1, typically.
2) Second-Order Partial Derivatives of x(ψ) w.r.t.φ: In this case, the last three terms in (60) vanish, again because of (61).We readily obtain, A 2-D coordinate system including the layout of the room where the measurements were taken is shown in Fig. 7. Also reported are the two selected positions p 1 and p 2 of (the center of gravity of) the Rx (virtual) array and the fixed position p Tx of the (single) Tx antenna.We recall that the mirror source method [74] computes the positions of predicted virtual sources associated with rays from the Tx antenna to the Rx array positions that undergo up to 5 reflections on walls or large objects (windows, boards).For the sake of conciseness we refer to virtual sources in the sequel as sources.The position, denoted by pl , of the predicted source corresponding to the lth SC, l = 1, . . .L detected by the algorithm is computed based on the estimated dispersion vector ψl using the relation pl = p + cτ l [cos( φl ) sin( φl )] T , where p either equals p 1 or p 2 .These positions are depicted in Fig. 7 as red crosses and blue pluses for the Rx array positions p 1 and p 2 , respectively.
The procedure described next attempts to associates detected sources and predicted sources.Possibly no, one, or more than one predicted SCs are associated to each detected SC as follows.A predicted SC is associated to a detected SC if their respective distances and angles are no more than, respectively, 10 cm (1/3 of the RRL in distance Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. in angle) apart.These selected values are within the same order of magnitude as, respectively, the 5 cm approximate accuracy of the floorplan (measured with a tape measure) and the CRLBs of the estimated distances and angles.The positions of successfully associated predicted sources are depicted in Fig. 7 as black circles and triangles for Rx positions p 1 and p 2 , respectively.
The algorithm is able to identify the LOS, most of the predicted first-order reflections and some predicted higher-order reflections for both Rx positions.Worth noting are the rays with reflections up to order five via the white board and the window highlighted in Fig. 7 and Fig. 5 and the second-order rays with reflections via the west plaster board and the east plaster board.The two former items are made of more reflective materials than the two latter.Furthermore, scattering from a metallic frame (see Fig. 5) that was not considered in the mirror source method could explain the detected source located at approximately [−8 6] m close to the west plaster board.For Rx position p 2 many detected sources are found in a region around [−15 8] m.They are likely to originate from scattering from this metallic frame.

Manuscript received 11
August 2023; revised 23 December 2023; accepted 18 February 2024.Date of publication 29 February 2024; date of current version 14 August 2024.This work was supported in part by the Austrian Science Fund (FWF) under Grant 10.55776/J4027; in part by the Christian Doppler Research Association; in part by the Austrian Federal Ministry for Digital and Economic Affairs; and in part by the National Foundation for Research, Technology, and Development.The associate editor coordinating the review of this article and approving it for publication was S. Mazuelas.(Stefan Grebien and Erik Leitinger contributed equally to this work.)(Corresponding author: Erik Leitinger.)

Fig. 1 .
Fig. 1.Layout of the array with its center of gravity p, the mth antenna element position p m and reference orientation o.The first SC originates from propagation along the direct path from the Tx to the Rx with angle φ1 and path distance cτ 1 .The kth SC is incident with angle φk and path distance cτ k .

8 end
their parameter estimates as well as the estimated parameters of the DC have converged.It then returns these converged values as the model estimates.

Fig. 4 .
Fig. 4. Detection and estimation of two closely spaced SCs in the scenario depicted in Subsubsec.VII-A4.Blue solid and red dashed lines refer to the bottom and top horizontal-axis, respectively.The panels in each row depict results obtained for the same SNR value, namely SNR = 10, 30, 50 dB, from the upper to the lower row.Panels in the columns depict from left to right the mean number of detected SCs, the relative frequency that exactly two SCs are detected, the RMSE of the estimated delays, the RMSE of the estimated angles, and the mean (unfilled) and RMSE (filled) of the moduli of the estimated amplitudes.

Fig. 5 .
Fig. 5. (a) Picture of the investigated room with the Rx, the Tx and some large-scale items labeled and (b) UWB antenna used at the Rx and the Tx.A floorplan of the room is given in Fig. 7 in Appendix D.
present the mean number of detected SCs ⟨ L⟩ and the relative frequency that exactly two SCs are detected ⟨1( L = 2)⟩ versus spacing in distance and angle, respectively.Both algorithms are able to reliably find the correct number of SCs provided the respective dispersion vectors of the two SCs are sufficiently apart.At high SNR (see Figs. 4(k)-(o)) the spacing values beyond which this occurs is as low as 0.15 m or 20 • for the system setting used in the study.Given that for this setting the RRL in distance is 0.3 m and that in angle is 56 • , this result demonstrates the superresolution capability of the proposed algorithm.At lower SNR (see Figs. 4(a)-(e)) these values rise towards the RRLs.Note that these values are not only influenced by AWGN but also by the DC, meaning that at high SNR the resolution

Fig. 6 .
Fig. 6.Selected application of the proposed algorithm to assess the dispersion characteristics (in delay and angle of arrival) of the channel from UWB SIMO measurement data: (a) estimated DAPS and dispersion vectors of detected SCs, (b) various estimated DPS, (c) residual DAPS and dispersion vectors of un-associated SCs, (d) various estimated APS. Detailed descriptions of the measurement setting and the depicted results are given in Subsec.VII-B.
cosine transfer function with roll-off factor 0.6 and bandwidth B = 1.6 GHz centered at f c = 6 GHz.The output signal with reduced bandwidth B is then Fourier transformed and sampled over [−B/2, +B/2] with frequency spacing ∆ = 6.8085MHz to produce a length N = 235 vector collecting these samples.A virtual 3 × 3 antenna array with 2 cm inter-element spacing is emulated by means of a single antenna mounted on a positioning table.Since the received signals are essentially noise-free, WGN was artificially added with power set such that SNR = 40 dB to emulate the model in (4), see Table 3.1] to 2|x(ψ • )| 2 and making use of [50, Theorem 4.4.1]combined with [50, Section 4.5.2]we obtain □ APPENDIX C COVARIANCE MATRIX OF THE χ 2 RANDOM FIELD 2|x(ψ • )| 2

3 )
Second-Order Partial Derivatives of x(ψ) w.r.t.τ and φ: In this case (61) make all terms in (60) vanish.Thus, E ∂x(ψ)∂x(ψ) OF THE SCS DETECTED BY THE ALGORITHM In this appendix, we provide a qualitative study that attempts to relate the SCs detected by the proposed algorithm to likely propagation mechanisms in the environment where the experimental data were collected.The results of this study supplement those presented in Subsec.VII-B.

Fig. 7 .
Fig.7.Floorplan of the investigated environment including the (fixed) Tx position p Tx , the two selected positions p 1 (red filled circle) and p 2 (blue filled triangle) of (the center of gravity of) the Rx array, the locations of estimated (virtual) sources for Rx array position p 1 (red crosses) and Rx array position p 2 (blue pluses) and the positions of associated predicted sources (black circles and triangles, respectively) computed with the mirror source model.