Dynamic Delay-Dispersive UWB-Radar Targets: Modeling and Estimation

This publication proposes a parametric data model and a gradient-based maximum likelihood estimator suitable for the description of delay-dispersive responses of multiple dynamic ultrawideband (UWB)-radar targets. The target responses are estimated jointly with the global target parameters range and velocity. The large relative bandwidth of UWB has consequences for model-based parameter estimation. On the one hand, the Doppler effect leads to a dispersive response in the Doppler spectrum and to a coupling of the target parameters that both need to be considered during modeling and estimation. On the other hand, the shape of an extended target results in a dispersive response in range, which can be resolved by the radar resolution. We consider this extended response as a parameter of interest, e.g., for the purpose of target recognition. Hence, we propose an efficient description and estimation of it by a finite impulse response (FIR) structure only imposing a restriction on the target’s dispersiveness in range. We evaluate the approach on simulations, compare it to state-of-the-art solutions, and provide a validation of the FIR model on measurements of a static scenario.

Abstract-This publication proposes a parametric data model and a gradient-based maximum likelihood estimator suitable for the description of delay-dispersive responses of multiple dynamic ultrawideband (UWB)-radar targets. The target responses are estimated jointly with the global target parameters range and velocity. The large relative bandwidth of UWB has consequences for model-based parameter estimation. On the one hand, the Doppler effect leads to a dispersive response in the Doppler spectrum and to a coupling of the target parameters that both need to be considered during modeling and estimation. On the other hand, the shape of an extended target results in a dispersive response in range, which can be resolved by the radar resolution. We consider this extended response as a parameter of interest, e.g., for the purpose of target recognition. Hence, we propose an efficient description and estimation of it by a finite impulse response (FIR) structure only imposing a restriction on the target's dispersiveness in range. We evaluate the approach on simulations, compare it to state-of-the-art solutions, and provide a validation of the FIR model on measurements of a static scenario.
Index Terms-Doppler effect, extended target model, highresolution parameter estimation, ultrawideband radar, velocity model.

I. INTRODUCTION
T HE range resolution of a radar system essentially depends on the bandwidth. A wide bandwidth is readily available at millimeter and subterahertz frequencies. Given a large bandwidth, the shape of an extended target shows up as a characteristic signature that is resolved in fast time (or range). However, some radar applications additionally require lower frequencies for manifold reasons. One advantage is the ability to penetrate materials, which improves with lower frequencies.
The frequency band of interest can range from a few hundred megahertz to several gigahertz where the latter is often termed ultrawideband (UWB). In this article, we assume a large relative bandwidth [2], [6] resulting from a large absolute bandwidth of several gigahertz at low frequencies.
For example, a relative bandwidth of more than 20% of the carrier frequency can be considered decisive for UWB.
Large relative bandwidths may cause issues with circuit design and algorithms. Considering signal parameter estimation, several algorithms are built on narrowband assumptions, i.e., model simplifications that are only valid for small relative bandwidths. A prominent example is the derivation of the velocity from the Doppler shift. Even with a wide bandwidth, one can still apply this principle of velocity estimation. However, this approach needs to be revised when the fractional bandwidth increases, as it is the case with UWB [7], [8]. Another example of invalid narrowband assumptions is that the target parameters range and velocity cannot be estimated separately anymore. Both effects are detailed in the sequel. In general, parameter estimation using wrong model assumptions becomes inherently biased since there is a mismatch between the model and the measured data [9], [10].
In traditional radar signal processing, the range-Doppler map, also known as scattering function [11], [12], [13], is estimated by a 2-D fast Fourier transform (FFT) assuming periodically repeated transmit signals. Typical examples are multicarrier sequences [orthogonal frequency-division multiplexing (OFDM)] or maximum length binary sequences (MLBSs). A matched filter response in the frequency domain transforms to fast time where the short correlation function indicates the target range. The respective Fourier transform along slow time yields the Doppler domain and the relative target speed is attributed to a Doppler shift resolved in delay and Doppler. However, from a physical point of view, the Doppler effect is a temporal scaling that can only be approximated as a frequency shift if the relative bandwidth of the waveform is small. Strictly speaking, the Doppler effect only results in a frequency shift for a sine wave.
Given UWB waveforms, a constant relative speed does not appear as a peak in the Doppler spectrum limited to a single range-Doppler resolution bin. It rather results in a dispersive response such that a Doppler shift-based model is not suitable for speed estimation. Hence, a modified formulation of the data model is required, which includes the relative velocity as 0018-926X © 2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
an explicit parameter. A further issue arises from the range migration of the impulse when the target is moving relative to the sensor. If the coherent integration time is chosen long enough (to increase the Doppler resolution), the pulse position migrates to the next range-Doppler resolution bin. Thus, the Doppler and range dimensions are no longer independent [14]. Finally, the Doppler effect additionally affects the observed characteristic signature of an extended target due to a time scaling of its response. If we want to model and evaluate this response, we need to know, on the one hand, the scaled waveform at the moment it hits the target. On the other hand, the target response to the scaled waveform is affected again by the Doppler effect on the way to the receiver. At the receiver, only the combined Doppler influence on the transmitting and receiving path is observable, which needs to be considered during target response estimation.
When it comes to modeling the characteristics of extended targets, current algorithms only allow this to a limited extent. Common multipath estimators (such as RIMAX [15]) assume a point-like interaction. Since an extended target spans a range of multiple resolution cells, it could be described by several reflection points. However, describing extended targets in this way would require an estimation of multiple and independent global range-velocity parameter pairs. Afterward, reflection points belonging to a single target have to be combined into a cluster of such pairs [16]. However, this approach has a drawback for the estimator. If the reflection points are closely spaced in the parameter domains, the estimated parameters are strongly correlated, resulting in high parameter variances for the individual reflection points. Alternatively, the model has to be adapted to directly incorporate extended target responses. Hence, we assume that a number of adjacent range samples (with the same velocity) belong to the same target. Furthermore, we assume that the response behaves linear and time-invariant (LTI) within the observation time. Consequently, we can select a finite impulse response (FIR) or infinite impulse response (IIR) approach as a structural (parametric) model of the target. If the extended target response is predominantly generated by independent scattering from the target's structure without internal multipath interaction, it can be modeled as an FIR structure. If the target's response is more strongly characterized by mutual interactions of scattering centers and resonances in cavities, IIR models may be more suitable.
Modeling the delay dispersion by such a compact, targetrelated model relaxes the overall estimation effort since one target is described by one global range-velocity parameter pair. Nevertheless, small variances in the target's trajectory and orientation will lead to a scintillation of the observed response resulting from reflections at structures larger than one wavelength. For point-like targets, this results in a fluctuation of amplitude and phase, requiring stochastic target descriptions such as Swerling models. Integration over time reduces the fluctuation and increases the probability of detection [17]. The response of an extended target will also be subject to a stochastic fluctuation due to the scintillation of individual scattering centers. Also, in this case, the fluctuations can be reduced by observing the target response over the coherent integration time. Furthermore, large-scale changes in the distance and orientation of the target can result in entirely different target responses. Hence, a physical interpretation of the model, e.g., for target recognition, requires the identification of poseinvariant features or matching to known angular-resolved target responses [18], [19].
In this article, we address two important topics: 1) relative target velocity estimation and 2) target response identification and modeling.
A. State-of-the-Art Many UWB radar applications have in common that either the sensing system or the radar targets are moving. An example for dynamic sensing systems is given in [20] where mobile and deployable sensing nodes cooperate in emergency situations, e.g., for environmental mapping and localization. Jovanoska et al. [21] and Zetik et al. [22] focused on dynamic targets, specifically on the detection and tracking of humans. Besides a variety of dynamic scenarios, some UWB applications do not only require a detection but also a target recognition [4], [20]. The contributions in [18] and [19] show that the extended target response in a highresolution range profile is valuable for target classification. However, extended target responses require new concepts for global target parameters as investigated in [23], e.g., for the global range of a human. To compensate Doppler spread and range migration in target velocity estimation, Dogaru [8] and He et al. [14] used matched filter and Keystone transformations for range-velocity map estimation. Nevertheless, a joint estimation of the global target parameters range and velocity as well as extended target responses is not yet analyzed.
Model-based target parameter estimation is investigated in a variety of contributions aimed at radio environments. The corresponding algorithms are RIMAX [15], [24] or the approaches presented in [25], [26], [27], and [28]. Although these algorithms provide a comprehensive description of the wave propagation by considering, e.g., angles of departure or arrival or bistatic antenna configurations, their applicability to UWB is limited by their narrowband assumptions. Relaxing these assumptions has as consequence that the target parameters cannot be estimated separately anymore since they are coupled, e.g., by the Doppler effect. Hence, Xu and Yarovoy [7] introduced coupling terms between range, velocity, and direction of arrival for their UWB target parameter estimator.
The aforementioned parameter estimators assume point-like targets. Hence, the underlying models do not account for any target-related spreads in the parameter domains such as range. Extended targets are therefore resolved into contributions of individual scattering centers. Building on these results, clustering algorithms are used to combine the contributions into clusters [29], [30], [31], [32]. Hence, a spread or dispersiveness in range is represented by a group of individual estimated scattering centers. However, prior estimation of individual scattering points can be a challenge for the estimator if the reflection points are closely spaced and, hence, differ very little in their parameters. Besides from model-based parameter estimation, Kuang et al. [33] investigated clustering algorithms on the radar detector's output to group extended target responses in the range-Doppler map.
Various mechanisms lead to a delay dispersion of a target's response. Some physical scattering events depend on frequency such as edge or corner diffraction. These mechanisms can be physically modeled by the geometrical theory of diffraction (GTD) [34], [35], [36], [37]. However, considering target parameter estimation in dynamic scenarios, this would require a comprehensive electromagnetic and angular resolved model, which does not generalize to different types of targets. Therefore, this approach is more suited for simulation studies.
If the target response is characterized by multipath reflections between scattering centers or strong resonant structures, it can be modeled as a superposition of damped exponential functions as in [38], [39], [40], and [41]. Here, the resonant structures are characterized by a resonance frequency and a damping factor describing the decay over time. Similarly, Saadane and Hayar [42] used an exponentially decaying pulse in time domain to account for resonant structures. Assuming that both, individual scattering centers and resonant structures, contribute to the target response, Eom and Chellappa [19] and Carriere and Moses [43] used autoregressive and moving average (ARMA) models known from time series analysis to model extended targets. This concept is similar to composing a target response of FIR (moving average) or IIR (autoregressive) structures. However, these approaches do not consider multitarget scenarios or global target parameters.
Delay-dispersive phenomena are also known for wideband radio channel modeling and estimation. A prevalent approach is to assign each propagation path an individual frequency transfer function (TF), which accounts for the dispersive properties of the path. Molisch [6] proposed to decompose the full TF into a number of subbands such that narrowband model simplifications are valid for each individual subband. This reduces the amount of necessary model parameters. Haneda and Takada [44] used a space-alternating generalized expectation-maximization (SAGE) algorithm to estimate individual parameters for each subband from channel measurements. In [45], a variant of the RIMAX algorithm is applied to a subband model, although parameters, such as time and direction of arrival, are kept constant over all bands. The advantage of these models is that they can capture a wide range of physical processes due to the generalized description of the TF. However, the approximations may hold for a narrow subband but still introduce discontinuities in the transition from one band to another, which are physically hard to justify. Hence, in [46], a sinc interpolation of smooth path TFs is introduced as an extension to the RIMAX algorithm.

B. Contributions
Building on the state-of-the-art, we propose a multitarget signal model (Section II) and a corresponding maximum likelihood (ML) parameter estimator (Section III) allowing a joint estimation of global target parameters (range and velocity) as well as the characteristic signature of a target. We model the delay-dispersive target responses as an FIR structure. This parametric representation also allows straightforward postprocessing of the target's response. Furthermore, we perform a comprehensive analysis of the approach by simulations (Section IV) and comparisons to state-of-theart solutions (Sections IV-B2, IV-B3, and IV-C1) and we verify the FIR approach on a static measurement example (Section V).
The proposed ML estimator is based on the RIMAX algorithm [15] and allows a grid-free estimation of the global target parameters. We discuss two major extensions to the RIMAX algorithm. First, a relaxation of the Kronecker-model assumptions is required since frequency and time are coupled by the Doppler effect. Second, the delay-dispersive model for each detected target needs to be estimated by a least-squares approach. Furthermore, the estimator includes heuristics based on statistical tests for an estimation of the number of present targets.

C. Notation
Let J x f denote the Jacobian matrix of a function f : R n → R m , which can also be a partial derivative, if f has more parameters than x. Let J 2 f denote the Hessian, i.e., the second-order derivatives, of the function f . We denote the Frobenius norm of a matrix A via ∥ A∥ F . Let 1 ∈ R n denote a vector of all ones, A H denote the Hermitian transpose of the matrix A, O denote the standard big-O notation, and z ∈ C denote the complex conjugate of z ∈ C. Let E denote the expectation operator.

II. DISPERSIVE MULTITARGET SIGNAL MODEL
We model the received signal R as response of a linear time-variant system to a periodically transmitted signal T This model requires that the received signal does not decorrelate with the transmitted signal due to the Doppler effect. This is valid, e.g., for short pulses or linear frequencymodulated sine waves [2, Ch. 2.2.6], [8], [47], though waveform design is not within the scope of this publication. Applying a matched filter to (II.1) provides an observation Y of the linear time-variant system S modeled by where N accounts for a yet-to-be-specified additive measurement noise process. The system itself is modeled as a superposition of P target responses S p The response of target p is composed of a TF as well as the global target parameters. These parameters consist of a delay τ p representing the traveled distance of the electromagnetic wave and a velocity v p . The velocity accounts for a timevariant delay with the speed of light c and is not equal to the absolute target velocity but the sum of its radial components relative Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
to the transmitter and receiver. In the case of a monostatic radar, the radial target velocity toward the radar is given by v p = v p /2. Considering the Fourier transform relation for a Dirac delta function we formulate the model for a single target response in complex baseband and time as where f c is the frequency employed for downmixing to baseband. The target's TF is represented by γ p : R → C. The three factors of (II.6) are further discussed in Sections II-A-II-C. The presented model separates the description of the target's TF via γ and τ from the velocity model parameterized by v. This involves that the target's TF is captured by the model as present in Y without considering the Doppler effect in the first place. We propose to compensate a Doppler-related scaling of γ after parameter estimation as derived in Section II-E and justified by numerical evaluations in Section IV-C2.

A. Target Frequency Transfer Function
We model the characteristic delay-dispersive signature of an extended target by an FIR structure. We represent this structure by a smooth estimate of the target's TF. The only prerequisite is that the signature has a limited width in time domain, allowing to describe the smooth frequency TF by a small number of sinc functions [46]. The sinc functions are centered at equidistantly spaced frequencies such that the model for a single target response γ p is formulated as where sinc(x) = sin(π x)/π x with sinc(0) = 1 and N s ∈ N.
In addition, we introduce a model-based bandwidth extension denoted by N e ∈ N 0 to mitigate boundary effects at the edges of the spectrum. This in analyzed in Section IV-B. Each sinc function is scaled by a constant complex weight γ ps . Considering the model for γ , one identifies two implications. First, we perform compression in the frequency domain by only estimating N s + 2N e amplitudes of the sinc functions instead of a dense frequency sampling. Second, the model provides a continuous and smooth representation of the TF. This process requires that γ is limited in the time domain to an interval of size T s . This limit is directly related by T s = 1 / S to the sinc function spacing S, which represents the frequency resolution of the model. This relationship is due to the fact that each sinc function is equivalent to a rectangular model window in the time domain with a width of T s . Due to this fact, a prior assumption about the maximum time extension of γ is required.

B. Simplification of the Delay Exponential
Due to the baseband representation of (II.6), the delay exponential consists of a linear phase accounting for the actual delay and a constant phase offset. This offset, however, is concentrated in time at τ = 0. Since we can reasonably assume that no target response is present at τ = 0, we can neglect this offset in the following.

C. Velocity and Doppler Model
From the velocity exponential in (II.6), we definê α : which is affine linear in f with proportionality constant v p and similar to the definition of the Doppler effect on a single frequency if v p ≪ c. Now, we slightly modify (II.8) by the following substitution [8]: This allows a significant interpretation. Whenever it holds that max( f ) ≪ f c , i.e., the signal covers a small bandwidth compared to the carrier frequency, the frequency dependence and, hence, Doppler scaling are negligible and the model is dominated by a common Doppler shift. If the bandwidth is not sufficiently small,α( f, α) accounts for a time-frequency coupling by introducing a time-dependent linear phase in the exponential in (II.6) besides the Doppler shift. This reveals that the model is most important for systems with high relative bandwidths, i.e., UWB systems. In the following, we consider α from (II.9) as a parameter implicitly estimating v by (II.10).

D. Discretized Model
To formulate a discrete observation model, we define a uniformly sampled version of S over frequency and time. We define these samples as where f > 0 and t > 0 are the sampling intervals in frequency and time, respectively, and N f , N t ∈ N are the numbers of acquired samples. The model parameters are continuous in the intervals τ p ∈ [0, 1/ f ) and α p ∈ [−0.5/ t, 0.5/ t). For the ease of implementation, we consider the normalization Summarizing, using (II.6), (II.7), and (II.10), our discrete model results in S ∈ C N f ×N t whose elements S i j are defined as where γ ∈ C P×N s +2·N e , τ ∈ R P , and α ∈ R P . Furthermore, we introduce the relative bandwidth

E. Doppler Effect on the Target's TF
The observed target TF γ over frequency (or fast time) will also be subject to the Doppler effect, which is not yet considered in (II.12). This model is only considering the modulation over (slow) time. Therefore, the TF is estimated in the scaled state and we propose to compensate it afterward as justified by simulations in Section IV-C2. To evaluate the influence on fast time, the Doppler effect can be modeled as a frequency mapping defined by considering the Doppler shift from the discretized version of (II.10). Hence, (II.14) describes where a frequency f i will be shifted due to the target's velocity represented by α. Since f i and f ′ i are normalized differently than α, η converts between the two normalizations.
The variable α ∈ [−0.5, 0.5] is normalized to the sampling interval in time t such that α = ±0.5 corresponds to the maximum unambiguous velocity ± v m . The Doppler shift in Hz corresponding to α is given by (II.9). Dividing it by f converts it to the normalization of The factor 2 assures that α = ±0.5 corresponds to the Doppler shift generated by v m . Using the definition of the relative bandwidth from Section II-D, f c can be expressed as Hence, the normalization conversion η reads as The Doppler scaling can be compensated by evaluating the TF-related parts of the model on a nonuniform frequency grid defined by (II.14) after the actual parameter estimation as To acquire the actual target response, f ′ i maps the frequency axis f i ∈ [−N f /2, N f /2 − 1] to the frequency values of the scaled target TF, i.e., to which frequency f i has been shifted due to the Doppler effect. In addition, the evaluation on the nonuniform frequency grid f ′ i leads to an offset of the delay parameter. To preserve the estimated global target delay, this offset is compensated by multiplying an additional delay exponential to the model in (II.18) with In Section III, we describe how we can estimate the model parameters γ , τ , and α from measured data by means of a ML estimator.

III. ML ESTIMATOR
The main computational obstacle of the proposed wideband data model is the coupling of frequency and slow time by the velocity model, which is made indispensable by large relative bandwidths. Since we assume to have noisy observations of the wave propagation in frequency and time as specified in (II.2), we need to consider the coupling in the estimation procedure of the parameters. Alternatively, one could integrate a model for the transmitted signal into the estimation routine allowing to exploit the wideband ambiguity function of the signal. In this way, the coupling could be realized by a 2-D matched filter for the transmitted signal directly transforming the received signal into the delay-velocity domain [8], [48]. This, however, would render the model and estimator waveform-specific.
For the presented approach, we draw inspiration from the RIMAX algorithm [49] and aim to provide an approach for extended target responses and velocities that can be viewed as an extension of it. The estimation problem can be split up in three intertwined subproblems. The first is to gain a reasonable estimate of the number of targets present in the observation. In addition, one has to find their parameters as demanded by our data model. Finally, we also have to account for possible measurement noise and stochastic components of the observation and have to estimate the corresponding distributions.
For the sake of clarity, we denote the original RIMAX algorithm and its applications as the narrowband case where, e.g., the velocity can be estimated from the Doppler shift. In contrast, the contributions of this publication target at the (ultra)wideband case.

A. Statistical Model for the Observations
To formulate a suitable statistical model for our proposed approach, we assume that an observation consists of the superposition of P ∈ N distinct target responses for which we can make use of the data model in (II.2) and its discretized version (II.12). Hence, an observation Y ∈ C N f ×N t follows: where N ∈ C N f ×N t is a random vector representing measurement noise at the sampled frequencies and time instances. In our case, we assume that the entries of N are independent and identically distributed (iid) according to a zero-mean circularly symmetric complex Gaussian distribution with noise variance σ 2 . Accordingly, each observation Y is a complex Gaussian random variable with mean S(γ , τ , α) and identity covariance scaled by σ 2 . Based on these assumptions, the objective is to estimate P ∈ N, (γ , τ , α) ∈ C N s +2N e ×P × R P × R P , and also σ > 0 from an observation Y.

B. ML Estimation
One popular approach for determining parameters from noisy observations is the principle of ML estimation. Given the statistical model in (III.1) following a Gaussian distribution, we can formulate the negative log likelihood function λ, given a certain observation Y, as: 2) In order to extract the ML estimate (γ ,τ ,α), one has to determine the global minimum of λ. A necessary condition for minimality is that the score function Jλ vanishes at (γ ,τ ,α) [50, Ch. 9.2]. The score function is the first-order derivative of λ defined as Jλ : C P×N s +2N e × R P × R P → C P×N s +2N e × R P × R P and is given by Also, at the optimum, the so-called Fisher information matrix (FIM) F, which is a positive semidefinite matrix, represents the expectation of the curvature of λ and is defined by Given the assumption that Y is following a Gaussian distribution, we can analytically calculate F via the Slepian-Bangs formula [51] given by Both quantities turn out to be vital components during the development of an algorithm that finds minimizers of the negative log likelihood as outlined in Section III-C2.
The optimization problem posed by ML is inherently hard to solve since the objective function is nonconvex, so a global minimizer exploiting the smoothness of λ will only converge to a local minimum, which is not necessarily a global optimum. Hence, we have to find a practically and computationally feasible approximation of an ML estimator.
The following Section III-C will first explain the approach taken to formulate such an estimator in the narrowband case. Then, we use this as a baseline to extend it to the considered wideband case in Section III-D.

C. Narrowband Iterative ML
For the arguably simpler case of narrow bandwidth, the ML estimator presented in [15] can make use of several algorithmic simplifications that result from the corresponding narrowband assumptions. These include specular reflections at point-like targets enabling propagation paths between transmitter and receiver as well as a narrowband Doppler model only considering a Doppler shift. Due to this assumption, the model for an observation may be simplified tõ with weightγ ∈ C P . One of the properties ofS is that the parametersτ andα are decoupled in the sense that they appear in distinct sampled and complex exponentials and hence influence distinct data dimensions, namely, frequency and time. This decoupling is realized algebraically by the use of the Kronecker or outer product to form the model for a single propagation path in (III.6). In the following, we will use· to mark quantities related to the narrowband case.
Similarly as before, we can define the negative log likelihood function based onS viã F . (III.7) Also, we can define the score functionλ and the FIMF similarly as before. The RIMAX algorithm works with the modelS and the respective likelihoodλ and consists of four main building blocks: a global search step [49,Ch. 5.1.5], where a single path is added to the model, a local search procedure [49,Ch. 5.2], where currently found parameters are refined iteratively by a second-order smooth optimizer, an estimator for the noise distribution [49,Ch. 6], and finally a method for removing unreliable paths from the model [49,Ch. 5.2.7.] based on their estimation variance.
As we will see next, the only step that needs substantial changes is the global search for a new path-in our case more precisely for a new target-and how to refine it after inclusion into the model. To motivate these changes, we briefly describe the approach of the algorithm in [49] that works with the narrowband models in (III.6) and (III.7).
On the highest level of abstraction, the RIMAX algorithm iteratively increases the model complexity, i.e., increases P until no more reliable paths can be added. Hence, the first problem one encounters is to find a suitable initial guess for the parameters of a newly added path. A good initial guess in step k is found, if the resulting parameter tuple (γ k ,τ k ,α k ) ∈ C k × R k ×R k is close enough to the optimum such that the objectivẽ λ is locally convex and the region of convexity contains the global minimum. In this case, a local optimizer initialized with (γ k ,τ k ,α k ) can recover the optimum given this initial guess.

1) Narrowband Path Search:
The approach presented in [49] essentially implements a grid-based and sequentially executed matched filter, which makes use of the Kronecker structure in (III.6). To describe this method, we define a f : (III.8) Next, we define two sets of grid points (III.9) Now, assume in step k that we already have a previous estimate (γ k−1 ,τ k−1 ,α k−1 ) at our disposal. Then, the so-called residual in step k − 1 can be defined viã To find a new tuple (τ ,α), we first find the index i f , which picks out the maximum element i f in the nonnegative vector Based on i f , we calculate the so-called beamformed residual Then, we find the index i t , which picks out the maximum entry in in order to finally derive the grid-based matched filter estimate via As a last step, we appendτ k andα k toτ k−1 andα k−1 , respectively, to formτ k andα k and find the path weightsγ k by means of least squares as given in (A.1). The methodology is summarized in Algorithm 1.
The main observation to make is that steps (III.11)-(III.13) are only valid since the parametersτ andα influence distinct data dimensions due to the narrowband assumption. Hence, the matched filter output can be computed efficiently with O(N 2 f · N t + N 2 t ) floating-point operations. In contrast, O(N 2 f · N 2 t ) operations would have been necessary, if we evaluated the grid for the two parameters jointly. Note that in [49, Ch. 5.1.2.], the algorithm is explained in much more generality such that it also applies to a multiple-input-multiple-output (MIMO)channel observation.
2) Narrowband Parameter Refinement: In order to mitigate the resolution limit imposed upon the global search procedure using G f and G t , the RIMAX algorithm executes SAGE [52] for all paths currently existing in the model each time a new path is added. Since SAGE only optimizes the parameters (γ ,τ ,α) for a single propagation path while keeping the others fixed, its computational cost is similar to that of Algorithm 1. In addition, Fessler and Hero [52] presented some convergence criteria, which makes SAGE a valid approach from an estimation perspective.
Once all required new paths are added by means of Algorithm 1 and SAGE, a second-order and joint optimization of the whole set of path parameters is carried out by which constitutes the so-called Fisher-scored gradient method [53] with step size ε, where the steepest descent direction is corrected according to the knowledge of the curvature at the current point of iteration. The iteration is denoted by ℓ. This algorithm has good convergence properties and comes at a moderate computational cost since no secondorder derivatives are necessary. At this point, we are in the place to formulate the algorithmic scheme adapted to our proposed data model. The important thing to notice is that for our purposes, it makes no difference if the iteration in (III.14) is executed based onS andλ or on S and λ using the respective score functions and FIMs. Furthermore, the computational complexity of S may be increased compared toS, but not as substantially as the complexity of the search for a new pathor more specifically a new target target-over the joint grid for τ and α. To retain the algorithmic advantages, the main idea of the wideband estimator is to use the initialization presented in Algorithm 1 with a slight modification.
To this end, we define the residual similar to (III.10) via Once we have established an initial estimate based onS, we replace the model by S while keeping the parameters (τ k−1 ,α k−1 ) fixed and determine γ k−1 by means of leastsquares now based on S. Then, we carry out the steps in (III.11)-(III.13) based on this wideband residual instead of the narrowband instanceR to obtain estimates for (τ k ,α k ). This is outlined in Algorithm 2. In this way, we have avoided using the wideband data model, where delay and Doppler are inherently coupled and have circumvented the joint 2-D path search. The parameter refinement with SAGE, however, is performed on the joint 2-D grid.
Once this initial estimate for (γ k , τ k , α k ) is established, we update it according to (γ k,ℓ+1 , τ k,ℓ+1 , α k,ℓ+1 ) Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. which is the wideband and hence more complex version of (III.14). However, we only need to evaluate the model S for the currently found k-paths, which results only in a mild increase in computational complexity for S compared toS. Both procedures are summarized in Fig. 1 to visualize the differences.

D. Computational Aspects
The iteration in (III.16) is usually paired with a stopping criterion that evaluates the relative changes of the likelihood function λ. Hence, in total, the necessary quantities that need to be computed efficiently are the likelihood λ and hence also its main building block S.
For the gradient iteration, we directly need access to Jλ, which also makes JS a quantity necessary to compute. Finally, to employ the correction by means of the inverse FIM, we see that by computing JS, we have already overcome the largest computational obstacle for F.

E. Model Order Selection and Noise Estimation
Although the data model of the proposed estimator is inherently different than the one presented in [49,Ch. 5.2.7], we have made straightforward modifications to the algorithm such that the methods for both presented therein are directly applicable to our case.
Section IV will study various effects and ramifications of the proposed data model and estimator.

IV. NUMERICAL EVALUATIONS
In this section, we analyze the proposed model and estimator from a numerical point of view and provide examples and visualizations. After introducing the simulation of the dispersive target responses, we discuss the implications of the wideband model. In addition, we compare our approach to state-of-the-art algorithms, i.e., RIMAX [49] in original form and a piecewise frequency-flat subband approach as proposed in [6] or [45]. Furthermore, we analyze a limitation of the approach and conclude this section by a discussion about two different aspects of the Doppler effect.
Whenever we visualize the time or Doppler domain, we compute an FFT utilizing a Hanning window over frequency or (slow) time. In addition, we average over the not-visualized data dimension.

A. Description of the Simulation
We simulate a multitarget response. Each individual target is modeled as extended in a way that it is represented by a cluster of specular reflections of closely spaced scattering centers. Each cluster is based on the narrowband model from (III.6) such that the full simulation model results in At this point, we consider a static scenario such that Y is constant over j except for the noise. P indicates the number of targets and Q indicates the cluster size, i.e., the number of closely spaced specular components composing a target. The complex amplitude is modeled as where a pq , ϕ pq ∼ U(0, 1) are drawn independently of a uniform distribution. The delays τ pq ∼ N (τ p , ϱ 2 p ) are drawn independently of a Gaussian distribution whose statistics depend on the target p. The variable τ p is the global target delay and ϱ p sets the width of the response. White Gaussian noise with noise power level σ 2 is added by n i j .
To add a dynamic component to the simulation, we extent (IV.1) by a velocity model in accordance with (II.12) The velocity parameter α varies between targets but is constant for all components composing a target. The static and dynamic simulation parameters are summarized in Table I.  the delay parameters τ p . The complete TF in Fig. 2(b) does not allow conclusions about the frequency characteristics of the individual targets though the multitarget nature of the wideband model and estimator captures the TF of each target separately, as visualized in Fig. 2(c) and (d). These plots reveal the different target characteristics, which can now be analyzed individually. Each TF is approximated by ten sinc functions spread over the spectrum. To enhance the model fit, we further exploit the fact that the sinc functions are not limited in the frequency domain but spread over the whole spectrum. For this, we define overlapping sinc functions specified by N e > 0 in (II.12). This allows to further enhance the model fit, especially at the edges of the spectrum. In this way, the model is not purely limited to the observed bandwidth up to a certain degree. The influence of N e on the estimation of the given TFs is exemplified in Fig. 3 where Fig. 3(a) shows the model fit given N e = 0 and Fig. 3(b) shows the model fit given N e = 2, which is also the setting for Fig. 2. A lack of overlapping sinc functions clearly influences the fit within the frequency band resulting in a ripple of the estimated TF.

B. Time-Invariant Analyses
Hence, for wideband target responses, the use of overlapping sinc functions is of great benefit.

1) Target Superposition:
The nature of the model requires a clear separation of the impulse responses (IRs) of multiple targets in the parameter domains, either in delay or velocity. If such a separation is not given, the model can only capture the superposition of the individual target responses and the estimator wrongly estimates a single target.
To analyze the transition from two targets wrongly estimated as a single one to a clear separation, we perform numerical experiments with the simulations from (IV.1) and Table I. The second target is varied in its τ parameter starting from an overlap with target 1 (τ 2 = 0.1) up to a clear separation (τ 2 = 0.15). For each τ 2 , we perform ten experiments with independent noise realizations and average the residuals accordingly. Fig. 4 shows the residual of the data model fit. The x-axis represents the varying τ 2 parameter and, hence, the distance τ between the targets. It is normalized to the model window T s . For ease of interpretation, ϱ = (ϱ 1 + ϱ 2 )/2 is added to τ since the target IRs have a certain width, and hence, the actual distance is accordingly smaller. The y-axis displays the global mean squared residual as well as the residuals between each target estimate and the corresponding simulated TF.
For small distances within the model window T s , both signals can be captured by a single target model. Therefore, the global residual is at the noise level though the residual of target 1 is quite large since it captures the superposition of both Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. targets. When the distance gets close to the model window, the width of the superimposed IR exceeds T s and cannot be captured by a single model component anymore. However, adding a second target to the model leads to unreliable parameter estimates and is hence omitted by the estimator. This is due to the fact that the model for the first target is not placed exactly at parameter τ 1 but at a different position where it fits best to the superimposed IR and maximizes the likelihood the most. Adding a second target will then lead to an overlap of the two model windows. After increasing the distance further, both target models can be placed in a way that each model describes an individual simulated target. Then, all three residuals drop to the noise level resulting in a good concordance between model and data.
The experiments reveal that the proposed model cannot distinguish between an individual target response or a superposition of several responses if they are not separated in the considered parameter domains. A decomposition into individual model components will fail, e.g., if several targets in a static scenario are distributed over the angular domains but share the same distance to the radar system. For separation, model and estimator would need to consider the angular domains.
2) Comparison to RIMAX: In Section III, we already discussed the similarities and differences between the RIMAX algorithm and the proposed wideband extensions. In this section, we want to investigate the claim of the introduction that such a point-like target model is indeed able to capture delay-dispersive response at the expense of resulting Fig. 4. Estimation residuals for a two-target scenario evaluated over a varying distance in delay between the targets. in a number of individual and unrelated model components. Therefore, we apply RIMAX to target 2 from the previous simulation [see Table I and (IV.1)]. The results are visualized in Fig. 5(a) and (b) and visually resemble the results from the wideband model in Fig. 2(a) and (d). As model order, ten paths are estimated whose delay parameters are highlighted in Fig. 5(a) by vertical lines. The total number of free model parameters for the RIMAX solution sums up to 20 (10-D τ and 10-D γ ) and is, hence, slightly higher as for the wideband model with 13 free parameters (1-D τ and (8 + 4)-dimensional γ ).
Nevertheless, the predefinition of the pulsewidth for the wideband model directly results in the estimation of a single target response. The RIMAX solution does not require such a predefinition. Thereby, the estimated paths do not reveal if they belong to single or multiple targets. A separation into individual and extended target responses would require a grouping or clustering of the estimated paths as a postprocessing step.
3) Comparison to a Constant Subband Model: Instead of modeling the dispersive target response by a superposition of specular components (Section IV-B2), one could approximate the target's frequency TF by decomposing the wideband spectrum into N s narrowbands, i.e., piecewise frequency-flat subbands. This is comparable to the approaches in [6], [44], and [45]. In our experiment, each of these subbands has its own complex amplitude γ ps as a model parameter, whereas τ is constant over all subbands. This is closely related to the proposed wideband model in (II.7), by replacing the sinc functions with rectangular functions as The TF is parameterized by γ ∈ C P×N s and the width of each subband is given by S > 0.  Table I and (IV.1)]. The frequency plot reveals the step-like approximation of the spectrum, which leads to a noticeable structure in the residual rising at each transition between two consecutive subbands.

C. Time-Variant Analyses
When it comes to time-variant scenarios, two aspects of the Doppler effect have to be considered. On the one hand, Doppler influences the signal modulation over (slow) time as already discussed in Section II-C. Given a high relative bandwidth, the modulation is caused by a scaling of the signal's spectrum. This effect is already incorporated into the proposed dispersive model in (II.12). On the other hand, the observed target TF γ over frequency (or fast time) will also be subject to the Doppler effect. The latter will be discussed in Section II-E.
The modulation over slow time is exemplary visualized in Fig. 6(a) as Doppler spectrum using simulations from (IV.3) and Table I. The different spreads of the two-target responses result from the velocity-dependent Doppler scaling. Since this scaling is incorporated into the model, both targets can be described by the scalar parameter α, which is directly related to a velocity by (II.9).
1) Comparison to RIMAX: In contrast to the wideband model, the RIMAX algorithm in the original form is built on narrowband assumptions and does only consider a Doppler shift assuming small relative bandwidths. In Section IV-B2, we discussed that a delay-dispersive target response can be captured by the RIMAX algorithm as a superposition of several individual model components. However, the mismatch of the Doppler model for larger relative bandwidths cannot be compensated by a superposition of multiple components, as visualized in Fig. 6(b). As in the example from Fig. 5(a), the algorithm spreads several components (in this case 20) over the Doppler domain. Nevertheless, the residual reveals that the model is not able to fit the simulated data sufficiently. Furthermore, while the decomposition into several delay components is physically justifiable by extended target geometries and individual scattering centers, the decomposition into distributed Doppler components is purely caused by model mismatch and prevents a derivation of the actual target velocity.
Furthermore, the impact of the model mismatch does not only depend on the target's velocity v but also on the velocity's relation to the sampling interval in time t. Since α is normalized to the sampling frequency 1/ t, small velocities can already lead to a significant spread in the Doppler spectrum.
2) Doppler Effect on the Target's TF: To evaluate the necessity of a Doppler effect compensation on the target's TF, we first need to include the Doppler scaling into the simulation introduced in Section IV-A. Similar to the proposed Doppler compensation in Section II-E, we can simulate Doppler scaled target responses by evaluating the delay exponentials of a target's TF in (IV.3) on a nonuniform frequency grid f ′ such that Variable f ′ i maps the frequency axis f i ∈ [−N f /2, N f /2 − 1] to the frequency values of the unmodulated target TF, i.e., to the frequencies from where the values at f i originate. Converting (IV.6) to f ′ i yields As with the Doppler compensation, we need to compensate for an artificial delay offset by with mean cluster delay τ p . Finally, to compare the results in this section for a wide range of simulated velocities, we assume that v m is adapted to each experiment such that v = 0.4· v m . Hence, α equals 0.2 for all experiments.
Evaluation: With the following experiments, we evaluate the necessity of a Doppler compensation of the estimated TF. Therefore, we performed two dynamic simulations, one as described in Section IV-A and another one with the modifications from above considering a Doppler scaling over fast time. We use the simulation parameters from Table I for target 2 and identical values drawn from the cluster statistics such that both TFs are identical if no Doppler scaling would be applied over fast time. We consider the previous simulation from Section IV-A as ground truth and the Doppler scaled simulation from this section as input to the parameter estimator. Hence, the estimator estimates the model parameters according to the scaled target TF. Afterward, we calculate the compensated TF as described in Section II-E.
For the evaluation, we define two mean-squared errors which captures the error between the ground-truth data and the estimated and uncompensated model given the Doppler simulation from (IV.7) as input and which captures the error between the ground-truth data and the compensated model from (II.18). The difference between the two errors E u − E c is visualized in Fig. 7 over two varying parameters: 1) the target velocity via v m by keeping α = 0.2 fixed and 2) the sampling interval f , i.e., the frequency resolution of the system. The latter is normalized to f c such that The varying resolution requires that the number of sinc functions N s is adapted accordingly to keep the size of the model window in the time domain fixed. N s is adapted by The difference E u − E c is drawn as colored areas on a decibel scale. As expected, the accuracy of the uncompensated model drops for rising target velocities and a rising frequency resolution. In other words, the uncompensated model performs worse if the Doppler shift per frequency is significant relative to the frequency resolution of the system. However, for most of the Fig. 7, the difference stays below 0.1 dB as highlighted by the red dashed line. Below this line, the gain of a Doppler compensation can be considered negligible. Two exemplary data points provide combinations of N f and v, which still do not require a Doppler compensation.
We conclude that the necessity of a Doppler compensation over fast time depends on the application and radar system. Hence, we do not incorporate Doppler-related terms into the model for the target's TF. This would apparently increase the complexity of the model and estimation routine. If necessary, the Doppler modulation can be compensated as postprocessing step, as described in Section II-E.

V. MEASUREMENTS
To conclude our investigations, we want to validate the applicability to an extended target on exemplary static short-range measurement data acquired with a m:explore M-sequence radar from Ilmsens GmbH operating in baseband from 0.1-6 GHz. The transmitter and receiver ports are connected to two DRH20E ridged horn antennas from RF SPIN, which are mounted on a stand in horizontal polarization resulting in a quasi-monostatic and single-polarized radar configuration. The signal-to-noise ratio (SNR) of this system effectively limits the usable frequency band to 1.5-5.5 GHz. The radar and signal parameters are summarized in Table II. The test object-a wooden plate with four metal corners as feet-is placed on Styrofoam in 0.8-cm distance to the antennas. Direct sources of reflections such as walls are covered with absorbers, though measurements have not been performed in an anechoic chamber. Furthermore, we acquire  Fig. 8(a) and (c)], the IR is dominated by a strong reflection originating from the surface of the plate. On the right, in Fig. 8(b) and (d), however, an extended echo is apparent dominated by the reflections of the two metal feet on the front.
To define the number of sinc functions for the model, we assume that the pulsewidth is directly related to the physical dimensions of the target. Therefore, we select a model window covering a range of 27 cm corresponding to 1.8 ns of free-space propagation. Given this, the number of sinc functions equals Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
Furthermore, we select N e = 2. Fig. 8(c) and (d) shows the resulting model window and the estimated model, respectively. Since the model window covers the extended target echo and the estimated model shows a good agreement with the measurement data, we conclude that N s = 8 is a reasonable choice for this target. Hence, selecting N s by the physical target dimensions is appropriate if no resonance effects occur at the target. In addition, the pulsewidth indeed not only depends on the maximum physical dimension of the target but also on the orientation relative to the radar.

VI. CONCLUSION
In this article, we propose a multitarget signal model and a corresponding ML parameter estimator targeting UWB-radar applications. Our approach allows a joint estimation of global target parameters (range and velocity) as well as the characteristic dispersive signature of a target. We formulate the model and estimator as extensions to the RIMAX algorithm relaxing its model simplifications designed for narrower bandwidths.
The main impact of the high relative bandwidth is that typical model simplifications concerning the Doppler effect are not valid anymore. This has three major consequences. First, the Doppler effect results in a dispersive response in the Doppler spectrum requiring an explicit model for the relative target velocity instead of a Doppler shift-based approach. Second, the global target parameters range and velocity require coupling terms and cannot be estimated independently of each other, as performed by RIMAX. Third, the Doppler effect results in a time scaling of the extended target signature, which needs to be considered during target response modeling. However, the aforementioned points are not equally important for all applications. While the velocity model and the parameter coupling are already necessary for small relative velocities, the Doppler scaling of the target response is only significant for high velocities. Therefore, we propose to compensate this scaling after the parameter estimation if required.
Regarding target response modeling, we assume that the response is composed of contributions from individual closely spaced scattering centers. Instead of estimating these scattering points individually, we propose to use an FIR structure. This, however, requires prior knowledge about the extent of the target response in fast time, i.e., the model order of the FIR structure. With this model, we can efficiently estimate extended target responses from noisy data. Nevertheless, this representation requires a clear separation of target responses in the considered parameter domains. Otherwise, the superposition of individual target responses is treated as belonging to a single target.
There are further extensions to the model worth analyzing. A model order selection based on statistical tests could render an assumption about the time extent of the target responses unnecessary. Furthermore, this would allow to use individual model orders for each target response. When it comes to resonant target structures, an extension by a recursive model is required. The estimator could be based on the Yule-Walker equations from autoregressive time series modeling. In addition, an increase of the parameter space, e.g., by direction of arrival estimation, could improve the separation of individual target responses and would add further valuable global target parameters. Furthermore, angular resolution can help to reject multipath reflections and clutter, which is not considered by the estimator right now. Alternatively, the parametric representation of the target response could be used for target recognition or to distinguish between desired target responses and undesired clutter.
In this publication, we only considered a dispersive response in range. However, one could also imagine that a target is extended in other domains such as in velocity, e.g., caused by moving rotor blades of a drone. It is conceivable to expand the FIR approach also to this domain. Nevertheless, we expect a drawback in performance between the FIR approach applied to multiple domains and an estimation of point scatterers in the range-velocity domain and a subsequent clustering approach. This needs to be investigated.

APPENDIX GENERALIZED LEAST SQUARES
Considering (III.2) and (III.7) together with (III.6) and (II.12), we realize that the parameters γ play a different role than τ and α.
Assume that we are given estimates τ 0 and α 0 for delay and Doppler, respectively. The minimizers of are linear functionsg : C N f ×N t → C P and g : C N f ×N t → C P×N s +2·N e with Y →g(Y) and Y → g(Y), respectively. Interestingly, in both cases, they can be designed to constitute the best linear unbiased estimator (BLUE). If both the data Y and the parameter γ were vectors, the mapping g would be realized by the so-called Moore-Penrose inverse. To avoid index-heavy notation, we skip the derivation of analytic expressions forg and g with the note that they can be derived analogously as the expressions of the Moore-Penrose matrices.