Super-Resolution Estimation of UWB Channels including the Diffuse Component — An SBL-Inspired Approach

In this paper, we present an iterative algorithm that detects and estimates the specular components and estimates the diffuse component of single-input—multiple-output (SIMO) ultra-wide-band (UWB) multipath channels. Specifically, the algorithm super-resolves the specular components in the delay– angle-of-arrival domain and estimates the parameters of a parametric model of the delay-angle power spectrum characterizing the diffuse component. Channel noise is also estimated. In essence, the algorithm solves the problem of estimating spectral lines (the specular components) in colored noise (generated by the diffuse component 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 specular component shall be retained or pruned. By relying to results from extreme-value analysis the threshold of this condition is suitably adapted to ensure a prescribed probability of detecting spurious specular components. Studies using synthetic and real channel measurement data demonstrate the virtues of the algorithm: it is able to still detect and accurately estimate specular components, 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 specular components than the actual ones. Finally, the algorithm is shown to outperform a state-of-the-art super-resolution channel estimator.


I. INTRODUCTION
Future wireless communication technologies will support a variety of services with high quality requirements, addressing performance metrics such as reliability, ultra-low 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 radio-based simultaneous localization and mapping (SLAM) [2]- [4].These examples emphasize the reliance of these technologies on extended, accurate channel state information.High-performance feasible parametric multiantenna 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 S. Grebien, E. Leitinger, and K. Witrisal are with the Laboratory of Signal Processing and Speech Communication, Graz University of Technology, Graz, Austria, and Christian Doppler Laboratory for Location-aware Electronic Systems (e-mail: {stefan.grebien,erik.leitinger,witrisal}@tugraz.at).B. Fleury is with the Institute of Telecommunications, Vienna University of Technology, Vienna, Austria, (bernard.fleury@tuwien.ac.at).S. Grebien and E. Leitinger have equally contributed as first authors.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 superresolution 1 tools to estimate their parameters.Expectationmaximization 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 socalled 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-aposteriori (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-power-exponential [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]- [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]- [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]- [38].However, some specificities of our underlying model -namely a two-dimensional dispersion domain and unknown colored noise -prevent a direct application of the method, see Subsection IV-5 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]- [46].These algorithms differ in their specific design criteria, such as (i) the chosen sparsity-inducing hierarchical prior model, e.g.gamma-Gaussian [40]- [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]- [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 the previously mentioned (iterative) SSR methods that apply grid refinement techniques [29], [35] can be viewed as particular instances of SSR methods with continuousparameter learning, which adapt their inherent restricted range of the dictionary parameter during the iterations. 3he 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], [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 dense component (DC) [9].The DC incorporates diffuse components 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 delayangle-of-arrival (angle) domain.The contributions of this paper are as follows: • We model the impact of the DC and AWGN as a colored noise, so that the problem becomes that of line spectral estimation [52] 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 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 algorithm with that of a state-of-the-art combined detection and estimation scheme that relies on the information criterion derived in [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.For the sake of simplicity we assume horizontal-only propagation. 5he array at the Rx has M elements located at p m ∈ R 2 , m ∈ {1, . . ., M } ≜ M, see Fig. 1.Its center of gravity is p = M −1 M m=1 p m and its orientation determined by the angle o as depicted in the figure .Signals are represented by means of their complex envelope with respect to a center frequency f c .Under the plane-wave assumption, the signal at the output of the mth antenna element reads In this expression s t; τ, φ, p m = e j2πfcg(φ,pm) s(t − (τ − g(φ, p m ))) (2) where s(t) is the transmitted signal with bandwidth B and g φ, p m = [cos(φ) sin(φ)](p m − 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 m relative to 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 m (t), m ∈ M are independent WGN with doublesided power spectral density N 0 /2.
We see from (1) that sufficient conditions for this identity to be accurate are that (a) the plane wave assumption holds over the Rx array aperture (determined by {p 1 , . . ., p m }) 6 , i.e. the array is located far away enough from the Tx and the objects in the environment that notably contribute to multipath propagation, 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.
In this study we assume that the delay-angle spread function h(τ, φ) is the sum of the superposition of a finite number, say K, of spectral lines representing SCs and a (spread) function ν(τ, φ) describing the DC, i.e.
Inserting the decomposition (3) in (1) yields 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 non-dispersive, 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.

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) 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 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 Q w = σ 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 v 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., S(f ; τ, φ, p (m) ) = e j2πfcg(φ,pm) 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], [56].It follows from Assumptions a) and b) that the covariance matrix of v factorizes as ) with ⊗ denoting the Kronecker product [9], [56].The first factor is the spatial correlation matrix with being the array response.The second factor is the delay correlation matrix with c) Experimental evidence shows that the DPS typically exhibits an exponentially decaying tail [9], [53], [57] and a smooth onset [53], [58].This behaviour 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], [59].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 Q takes the following form: where η = [σ 2 P ϑ] with ϑ defined above and Q f given in (13) with p(τ ) = p(τ ; ϑ) according to (14).Hence, because of Assumption c),

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 {S( ψ)} ψ∈Ψ 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 with s(ψ l ) defined similarly to (6).Under the assumptions made in Subsection II, the likelihood function of this model reads (18) 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

B. Inference Method
The proposed method is inspired from SBL [20].First it computes a MAP estimate of ψ, η, and γ from the joint posterior of these random vectors; then it uses this 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 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 (23), and thus in ( 24) and (25).
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 ( 26) 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 (26) 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 (29) 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 ) 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 (30) contains a condition that determines when the lth spectral line shall be discarded (γ l = ∞).
4) Estimation of the Weights: Inserting the estimates ψ, η, and γ in ( 23) yields the Gaussian PDF with mean (see (24)) and covariance matrix (see (25)) that is used as an approximation of the posterior PDF of α.
In (38), 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.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], [60], we increase the initial threshold κ = 1 in the pruning condition (30) 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 (49).

Discussion
The updating step (27) 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 (27) 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 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, S(ψ) 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.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 parametric class 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-dimensional [36], only an approximate such program could be formulated to date for higher dimensional 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 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.
It is shown in [51] that as a result of the first part in the assumption the matrix Q v in ( 11) is centro-hermitian 13 [63] and therefore Q 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 s( ψl ) is nearly orthogonal to any columns of S( ψl) for each l = K + 1, . . ., L.
As a result of Assumption 2, as M N grows large, with high probability (31) and ( 32) 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 (30), with threshold set to κ is close to when M N is sufficiently large.Let us consider the circularlysymmetric complex Gaussian random field on Ψ defined as with ψ • ∈ Ψ.Then, (39) can be recast as Theorem 1.Under Assumption 1 we have the asymptotic equivalence Furthermore, Here, Λ(ψ • ) is the non-negative definite matrix given in (45), incorporates the impact of the array aperture, while b(τ ) and a(τ ) incorporate the impact of the frequency aperture (spectrum S(f )) and colored noise.
Proof.As shown in [51] it follows from Assumption 1 that the real and imaginary parts of the Gaussian field x(ψ • ) in (40) exhibits the following properties: a random field χ 2 on Ψ with two degrees of freedom [49], [50]. 14Note 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 14 The real and imaginary parts of √ 2 x(ψ•) have unit variance, in accordance with the definition of a χ 2 process.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.
The function of κ in the asymptotic equivalence (42) provides a tight approximation of P f (κ) 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 f (κ ⋆ ) close to ϵ.The next lemma essentially gives this inverse function.
We conclude from ( 41), ( 42) and ( 49) that P f (κ ⋆ (ϵ)) ≈ ϵ for ϵ sufficiently small and M N sufficiently large.Thus, the function κ ⋆ (ϵ) provides a means to control the probability of detecting spurious components provided M N is sufficiently large.We can use the following asymptotic behavior of the function W −1 (u) as u → 0 to obtain a tight approximation of κ ⋆ (ϵ) for ϵ small: Algorithm 1: Main Input : Measurement vector y Output: L, ψ, η, μ, and with σ2 ← ∥y∥ 2 N M , P ← 0, θ ← 0 (AWGN only), and initFlag ← true ζ l using ( 31) and ( 32) append ψl to ψ and γl to γ where o(•) denotes the little-o notation [65].Making use of this identity, the equality in (49) can be recast as The right-hand expression with the term o(ϵ) dropped provides a tight approximation of κ ⋆ (ϵ) when ϵ is sufficiently small, provided M N is sufficiently large.

A. Examples:
We illustrate the right-hand expression in (43) 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 (43) becomes 4π N 2 −1 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 candidates 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 their parameter estimates as well as the estimated parameters of the DC are converged.It then returns these converged values as the model estimates.
The initial iterations in the do-while-loop of Algorithm 1 are executed while considering measurement noise only, i.e. by using Q in (15) with P set to zero whenever Q 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 the estimates of the parameters of the DC DPS are updated, i.e. the full covariance matrix Q in ( 15) is accounted for in the update equations of Procedures 1 and 2.

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 Q 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 specular-to-dense-ratio SDR = 10 log 10 1 M ∥ k∈K αk s( ψ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 s( ψk )∥ 2 + P B)/σ 2 .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 used settings are as follows: N = 27, i.e. ∆ = 59.26MHz, and the array has dimension 3 × 3 with 2 cm inter-element spacing, i.e.M = 9.Our algorithm uses a fixed threshold κ * = 4 and L = 50.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 Q 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 (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 (42) 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 parameter 2η = [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 ( 41) and ( 52) as approximations of the probabilities of detecting spurious SCs and missing an active SC, respectively, we modify the settings of the simulation scenario in Subsection VII-A1 as follows: The array has size 5 × 2; N = 54; L = 10; SNR = {5, 10, 20} dB; 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 f (κ) in (41) (dashed lines) and P m (κ) in (52) (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 [66].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 Q in (10) used in the generative model: Q has the simplified form (15) (Fig. 2a and Fig. 2b) and Q has the general form (10) (Fig. 2c). 16Furthermore, under the first assumption we distinguish between the two cases where the matrix Q 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 Q used 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 in Subsection VII-A2 as follows: the array is linear, i.e. 1 × M , with inter-element spacing of 2 cm; N = 27; L = 50; SNR = 20 dB; P = 0 (WGN only); ϵ = 10 −2 ; σ 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. 17

4) 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 vector.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 • ) [67].The complex amplitudes of both SCs have unit magnitude and their respective phases are drawn uniformly and independently.Other system param-eters are set as follows: the array dimension is 3 × 3 with 2 cm inter-element spacing; SDR = 6 dB; SNR = {10, 30, 50} dB; N = 54; ϵ = 10 −3 ; L = 10.
In Fig. 4, we compare the performance of the proposed algorithm (circles) with a maximum-likelihood estimation (MLE) algorithm (triangles) inspired by [9].The modification consists in adapting the algorithm to our application scenario, adopting the same scheduling and the same thresholding as in our algorithm.Apart from the used scheduling, the modified algorithm strongly resembles that in [14] with the information criterion adapted to our scenario.Performance versus spacings ∆d and ∆φ are depicted as, respectively, blue solid and red dashed curves.
The first two columns of Fig. 4 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 Fig. 4 (k) to (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 Fig. 4 (a) to (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 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 [68].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 meansquare error (MMSE) estimate based on the hyperparameter estimates, see (37) and (38).When the estimated dispersion vectors of the two detected SCs are closely spaced, the leastsquares 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 [69].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-raisedcosine transfer function with roll-off factor 0. The room where the measurements were performed is depicted in Fig. 5. Based on a layout of it, see Fig. 7 in [64, Appendix B], a classical mirror source method [70] 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 details in [64,Appendix B].
The following analysis concerns measurements obtained with Rx position p 1 depicted in Fig. 7 in [64, Appendix B].Fig. 6(a) depicts the estimated DAPS computed from the received signal [70].Note that this power spectrum incorporates the smoothing function of the aperture of the measurement equipment [62].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

VIII. CONCLUSIONS
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 diffuse component plus noise in an ultra-wide band single-inputmultiple-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 diffuse 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 his high efficiency the proposed algorithm has promising potential applications in all aspects of wireless communications that exploit extended channel state information, such as integrated sensing and communications (ISAC) and radio-based localization.
We now proceed with the computation of the entries of (45).
1) Second-order Partial Derivatives of x(ψ) w.r.t.τ : One can easily check that because of (54) the second and third terms of (53) vanish in this case.As a result, we can write and is a delay-dependent loss factor that depends of the structure of the noise vector n.Note that if n is white, e(τ ) = 1, otherwise e(τ ) < 1, typically.
2) Second-order Partial Derivatives of x(ψ) w.r.t.φ: In this case, the last three terms in (53) vanish, again because of (54).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 [70] 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 of the SC 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) and 5 • (1/10 of the RRL 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 higherorder 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 secondorder 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, that are likely to originate from scattering from this metallic frame.

Fig. 1 .
Fig.1.Layout of the array with its center of gravity p, the mth element position pm and reference orientation o.The 1th 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 .

Fig. 2 .
Fig.2.Comparison of the probalilities P f (κ) and Pm(κ) 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.

Fig. 3 .
Wideband versus narrowband detection and estimation of SCs in the AWGN channel: (a) mean number of detected SCs, (b) 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 .

Fig. 4 .Fig. 5 .
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.
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).The algorithm uses the following settings: L = 50; ϵ = 10 −2 .

Fig. 6 .
5 • (1/10 of the RRL), see [64, Appendix B.A] 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 Selected application of the 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.computed from the residual signal y−S( Ψ) μ.It also includes the dispersion vectors of detected SCs (blue triangles) 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.

3 )
Second-order Partial Derivatives of x(ψ) w.r.t.τ and φ: In this case (54) make all terms in (53) vanish.Thus,E ∂x(ψ)∂x(ψ) * ∂τ ∂φ = 0.(59)APPENDIX B VALIDATION OF THE SMCS DETECTED BY THEALGORITHMIn this appendix, we provide a qualitative study that attempt to relate the SCs detected by the proposed algorithm to probable 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).