Joint Transmit and Receive Filter Optimization for Sub-Nyquist Delay-Doppler Estimation

TL;DR

This paper introduces a joint transmit and receive filter optimization framework for delay-Doppler estimation that exceeds traditional bandwidth constraints, improving parameter estimation accuracy through Bayesian Cramér-Rao bound optimization.

Contribution

It proposes a novel transceiver design that violates the Nyquist-Shannon bandwidth limit by jointly optimizing filters for better delay-Doppler estimation accuracy.

Findings

Optimized filters outperform traditional designs in estimation accuracy.

The approach explores the Pareto-optimal filter region under weighted MSE.

Monte-Carlo simulations confirm improved estimator performance.

Abstract

In this article, a framework is presented for the joint optimization of the analog transmit and receive filter with respect to a parameter estimation problem. At the receiver, conventional signal processing systems restrict the two-sided bandwidth of the analog pre-filter $B$ to the rate of the analog-to-digital converter $f_s$ to comply with the well-known Nyquist-Shannon sampling theorem. In contrast, here we consider a transceiver that by design violates the common paradigm $B\leq f_s$. To this end, at the receiver, we allow for a higher pre-filter bandwidth $B>f_s$ and study the achievable parameter estimation accuracy under a fixed sampling rate when the transmit and receive filter are jointly optimized with respect to the Bayesian Cram\'{e}r-Rao lower bound. For the case of delay-Doppler estimation, we propose to approximate the required Fisher information matrix and solve the transceiver design problem by an alternating optimization algorithm. The presented approach allows us to explore the Pareto-optimal region spanned by transmit and receive filters which are favorable under a weighted mean squared error criterion. We also discuss the computational complexity of the obtained transceiver design by visualizing the resulting ambiguity function. Finally, we verify the performance of the optimized designs by Monte-Carlo simulations of a likelihood-based estimator.

Figures

Equations

$$\breve{x}(t)$$

$$\breve{y}(t)=\mathcal{R}\{\breve{x},\mathbf{\theta }\}(t)+\breve{\eta }(t).$$

$$y(t)$$

$$=v(t;\mathbf{\theta })+\eta (t).$$

$$\mathbf{y}=\mathbf{v}(\mathbf{\theta })+\mathbf{\eta },$$

$$\mathbf{u}$$

$${\mathbf{R}}_{\mathbf{\eta }}={\mathrm{E}}_{\mathbf{\eta }}[\mathbf{\eta }{\mathbf{\eta }}^H]\in {\mathbb{C}}^{N\times N}.$$

$$\mathrm{MSE}(\boldsymbol{M})=\mathrm{tr}\big{(}\boldsymbol{M}\boldsymbol{R}_{\boldsymbol{\epsilon}}\big{)},$$

$$\boldsymbol{R}_{\boldsymbol{\epsilon}}=\mathrm{E}_{\boldsymbol{y},\boldsymbol{\theta}}\Big{[}\big{(}\boldsymbol{\hat{\theta}}(\boldsymbol{y})-\boldsymbol{\theta}\big{)}\big{(}\boldsymbol{\hat{\theta}}(\boldsymbol{y})-\boldsymbol{\theta}\big{)}^{\operatorname{\operatorname{T}}}\Big{]}.$$

$${\mathbf{R}}_{\mathbf{ϵ}}⪰{\mathbf{J}}_B^{-1},$$

$${\mathbf{J}}_B={\mathbf{J}}_D+{\mathbf{J}}_P.$$

$$\boldsymbol{J}_{D}=\mathrm{E}_{\boldsymbol{\theta}}\big{[}\boldsymbol{J}_{F}(\boldsymbol{\theta})\big{]},$$

$$[\boldsymbol{J}_{F}(\boldsymbol{\theta})]_{ij}=-\mathrm{E}_{\boldsymbol{y}|\boldsymbol{\theta}}\Bigg{[}\frac{\partial^{2}\ln p(\boldsymbol{y}|\boldsymbol{\theta})}{\partial[\boldsymbol{\theta}]_{i}\partial[\boldsymbol{\theta}]_{j}}\Bigg{]}.$$

$$[{\mathbf{J}}_F(\mathbf{\theta }){]}_{ij}=2\mathrm{Re}\left\{{\left(\frac{\partial \mathbf{v}(\mathbf{\theta })}{\partial [\mathbf{\theta }{]}_i}\right)}^{\mathrm{H}}{\mathbf{R}}_{\mathbf{\eta }}^{-1}\left(\frac{\partial \mathbf{v}(\mathbf{\theta })}{\partial [\mathbf{\theta }{]}_j}\right)\right\}.$$

$$[\boldsymbol{J}_{P}]_{ij}=-\mathrm{E}_{\boldsymbol{\theta}}\Bigg{[}\frac{\partial^{2}\ln p(\boldsymbol{\theta})}{\partial[\boldsymbol{\theta}]_{i}\partial[\boldsymbol{\theta}]_{j}}\Bigg{]}.$$

$$\underset{g(t),h(t)}{\mathrm{argmin}}\mathrm{tr}(\mathbf{M}{\mathbf{R}}_{\mathbf{ϵ}}),s.t.\frac{1}{T_0}{\int }_{T_0}|\breve{x}(t){|}^2\mathrm{d}t\le P_T.$$

$$\underset{g(t),h(t)}{\mathrm{argmin}}\mathrm{tr}(\mathbf{M}{\mathbf{J}}_B^{-1}),s.t.\frac{1}{T_0}{\int }_{T_0}|\breve{x}(t){|}^2\mathrm{d}t\le P_T.$$

$$\underset{g(t),h(t)}{\mathrm{argmax}}\mathrm{tr}({\mathbf{M}}^{\prime }{\mathbf{J}}_B),s.t.\frac{1}{T_0}{\int }_{T_0}|\breve{x}(t){|}^2\mathrm{d}t\le P_T,$$

$$\underset{g(t),h(t)}{\mathrm{argmax}}\mathrm{tr}({\mathbf{M}}^{\prime }{\mathbf{J}}_D),s.t.\frac{1}{T_0}{\int }_{T_0}|\breve{x}(t){|}^2\mathrm{d}t\le P_T.$$

$$\mathcal{R}\{\breve{x},\mathbf{\theta }\}(t)=\gamma \breve{x}(t-\tau ){\mathrm{e}}^{\mathrm{j}2\pi \nu t},$$

$$\mathbf{\theta }={\left[\begin{array}{cc}\tau & \nu \end{array}\right]}^{\mathrm{T}}.$$

$$\mathbf{\theta }\sim \mathcal{N}(0,{\mathbf{R}}_{\mathbf{\theta }})$$

$${\mathbf{R}}_{\mathbf{\theta }}=\left[\begin{array}{cc}{\sigma }_{\tau }^2& 0\\ 0& {\sigma }_{\nu }^2\end{array}\right].$$

$$\breve{x}(t)=\sum _{k=-\frac{K}{2}}^{\frac{K}{2}-1}{G}_k{\mathrm{e}}^{\mathrm{j}k{\omega }_{0}t},$$

$$G_k=\frac{1}{T_0}\underset{-\infty }{\overset{\infty }{\int }}g(t){\mathrm{e}}^{\mathrm{j}k{\omega }_{0}t}\mathrm{d}t.$$

$$\operatorname{\mathcal{K}}=\Big{\{}-\frac{K}{2},-\frac{K}{2}+1,\dots,\frac{K}{2}-1\Big{\}}$$

$$v(t;\mathbf{\theta })$$

$$=\gamma {\mathrm{e}}^{\mathrm{j}2\pi \nu t}\sum _{k\in K}{\mathrm{e}}^{\mathrm{j}k{\omega }_{0}t}{\mathrm{e}}^{-\mathrm{j}k{\omega }_0\tau }{G}_{k}H(k{\omega }_0+2\pi \nu ),$$

$$H(\omega )=\underset{-\infty }{\overset{\infty }{\int }}h(t){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

$$v(n{T}_s;\mathbf{\theta })=\gamma {\mathrm{e}}^{\mathrm{j}2\pi \nu n{T}_s}\sum _{k\in K}{\mathrm{e}}^{\mathrm{j}2\pi \frac{kn}{N}}{\mathrm{e}}^{-\mathrm{j}k{\omega }_0\tau }{G}_{k}H(k{\omega }_0+2\pi \nu ).$$

$${\mathcal{F}}_k=\left\{k-N{L}_k^{-},\dots ,k,k+N,\dots ,k+N{L}_k^{+}\right\}.$$

Full text

I Introduction

The inference of unknown parameters is of interest in technical applications such as radar, sonar, image analysis, biomedicine or seismology. In radar systems knowledge of the delay-Doppler shift can be used to determine the distance and velocity of a target object. For imaging and biomedical applications conclusions about the test medium can be drawn from the received signal parameters. During the design of such systems, resource and hardware restrictions, e.g., power consumption and computational effort, have to be taken into consideration and a favorable trade-off between high system performance and low system complexity has to be found.

At the transmitter, using large bandwidths is, in general, favorable for the performance. Concerning the analog hardware costs at the receiver side, in particular, the sampling rate of the analog-to-digital (A/D) conversion has been identified as a bottleneck [1, 2]. In general, one can take two different engineering perspectives regarding the performance-complexity trade-off associated with the digitization process at the receiver. On the one hand, if the sampling rate $f_{s}$ is limited, fulfilling the sampling theorem by setting the analog pre-filter bandwidth to $B\leq f_{s}$ means that the performance is bounded by the complexity of the sampling device. On the other hand, if a certain bandwidth $B$ is available to the transmitter, fulfilling the Nyquist-Shannon theorem requires one to provide sufficiently large power and hardware resources at the receiver to perform sampling with $f_{s}\geq B$. While choosing $B\leq f_{s}$ guarantees perfect reconstruction of the band-limited analog received signal from its digital samples, for the signal processing task of pilot-based estimation we are instead interested in the accurate recovery of only a few signal parameters using a probabilistic approach. Since high estimation accuracy rather than a low reconstruction error is the desired goal, compliance with the sampling theorem can be relaxed. This leads to the fundamental question of how to design signal parameter estimation methods and systems when commonly used principles like the Nyquist-Shannon theorem are set aside.

II Problem Formulation

Before treating the special case of delay-Doppler estimation, we state the considered technical problem generically.

II-A Generic System Model

Consider the system model in Fig. 1,

where a periodic analog waveform $\breve{x}(t)\in\mathbb{C},t\in\mathbb{R}$ is emitted by continuously repeating the waveform $g(t)\in\mathbb{C}$ of bandwidth $B$,

[TABLE]

in intervals of $T_{0}$. The transmit signal (1) is transformed according to a parametric propagation model $\mathcal{R}\{\breve{x},\boldsymbol{\theta}\}(t)$, which is characterized by a functional operator with parameters $\boldsymbol{\theta}\in\mathbb{R}^{D}$. At the receiver, the waveform is perturbed by additive white Gaussian noise $\breve{\eta}(t)\in\mathbb{C}$ with constant power spectral density $N_{0}$, resulting in the analog sensor signal $\breve{y}(t)\in\mathbb{C}$,

[TABLE]

The signal (2) is filtered by a linear, time-invariant filter $h(t)\in\mathbb{C}$, such that the final analog received signal is given by

[TABLE]

Note that here the filtering operation with $h(t)$ is performed in the analog domain before the sub-Nyquist A/D conversion. Therefore, it is possible to shape the aliasing effect and diminish the associated information loss. In contrast, a digital receive filter would act on the sampled data and consequently could not recover the loss induced by sampling artifacts.

The signal $y(t)\in\mathbb{C}$ is sampled in intervals of $T_{s}=\frac{1}{f_{s}}$ for the duration $T_{0}$, resulting in $N=\frac{T_{0}}{T_{s}}\in\mathbb{N}$ samples

[TABLE]

where we use vector notation to denote $N$ samples

[TABLE]

of an analog waveform $u(t)\in\mathbb{C}$ with $t\in[-\frac{T_{0}}{2},\frac{T_{0}}{2})$. We assume that the noise samples $\boldsymbol{\eta}\in\mathbb{C}^{N}$ in (4) follow a multivariate zero-mean Gaussian distribution with covariance

[TABLE]

II-B Parameter Estimation Task

Throughout this work a pilot-based Bayesian estimation approach is assumed, i.e., the receiver has perfect knowledge about the transmit signal $\breve{x}(t)$, the receive filter $h(t)$, the sampling rate $f_{s}$ and the functional form of $\mathcal{R}\{\breve{x},\boldsymbol{\theta}\}(t)$ for each $\boldsymbol{\theta}$. Therefore, the exact probabilistic model $p(\boldsymbol{y}|\boldsymbol{\theta})$ is available at the receiver, while the propagation parameters $\boldsymbol{\theta}$ are assumed to be unknown random variables distributed according to a known probability law $\boldsymbol{\theta}\sim p(\boldsymbol{\theta})$.

The goal of the receiver is to infer the unknown parameters $\boldsymbol{\theta}$ based on the received digital data $\boldsymbol{y}$ using an appropriate estimation algorithm $\boldsymbol{\hat{\theta}}(\boldsymbol{y})$. The weighted estimation error of the processing procedure $\boldsymbol{\hat{\theta}}(\boldsymbol{y})$ is

[TABLE]

for some positive semidefinite weighting matrix $\boldsymbol{M}\in\mathbb{R}^{D\times D}$, where we define the mean squared error (MSE) matrix by

[TABLE]

As a direct computation of the MSE matrix (8) is generally intractable, a performance analysis is usually based on the Bayesian Cramér-Rao lower bound (BCRLB) [3, p. 5]

[TABLE]

where $\boldsymbol{J}_{B}$ is the Bayesian information matrix (BIM)

[TABLE]

The BCRLB (9) forms a fundamental limit for the achievable estimation accuracy (8). The first term on the right-hand side of (10) represents the expected Fisher information matrix (EFIM)

[TABLE]

with the Fisher information matrix (FIM) exhibiting entries

[TABLE]

The matrix (11) characterizes the average information about the parameters $\boldsymbol{\theta}$ that is embodied in the data $\boldsymbol{y}$. For the Gaussian model (4), the FIM entries (12) are given by [4, p. 525]

[TABLE]

The second summand in (9) is the prior information matrix (PIM) with entries

[TABLE]

It specifies the information about the unknown parameters which is contained in the prior $p(\boldsymbol{\theta})$.

Accessing and optimizing the system performance (8) based on a theoretical measure like (10) has the advantage that the achievable accuracy level can be computed analytically without extensive simulations of the algorithm $\boldsymbol{\hat{\theta}}(\boldsymbol{y})$.

II-C Transceiver Design Problem

The design problem of finding analog transmit and receive filters $(g^{\star}(t),h^{\star}(t))$ that minimize the MSE (7) of the estimation algorithm $\boldsymbol{\hat{\theta}}(\boldsymbol{y})$ under a particular weighting $\boldsymbol{M}$, subject to a transmit power constraint $P_{T}$, can be formulated as

[TABLE]

For simplicity, we assume that the system operates in the linear regime of the transmit power amplifier and therefore do not include a peak power constraint.

Although the BRCLB (9) can be achieved with equality only under special conditions [3, p. 5], it closely characterizes the trend of the estimation performance. This is illustrated by simulation of the MSE in Section XI. It is therefore possible to formulate (15) based on the BIM (10)

[TABLE]

To avoid optimization with respect to $\boldsymbol{J}_{B}^{-1}$ in (16), we consider an alternative maximization problem

[TABLE]

which is equivalent to (16) in a boundary preserving sense. It can be shown that if $(g^{\star}(t),h^{\star}(t))$ is a solution of the maximization problem (17) with $\boldsymbol{M}^{\prime}$, there exists a weighting matrix $\boldsymbol{M}$ (not necessarily equal to $\boldsymbol{M}^{\prime}$) for which the original optimization problem (16) has the same solution $(g^{\star}(t),h^{\star}(t))$ [5]. Since $\boldsymbol{J}_{P}$ is independent of $g(t)$ and $h(t)$, (17) simplifies to

[TABLE]

III Contribution and Outline

For the example of a single-input single-output (SISO) system with delay-Doppler shift, we develop a framework to solve the transceiver optimization problem (18) under the sub-Nyquist condition $B>f_{s}$. Note that such pilot-based delay-Doppler estimation problems are key for wireless applications like radar or radio-based localization and positioning.

After providing an overview on related work (Section IV), we

• •

Introduce the receive model with delay-Doppler shift and derive its FIM under arbitrary transmit and receive filters (Section V).

• •

Show how to approximate the information matrix $\boldsymbol{J}_{F}(\boldsymbol{\theta})$ and $\boldsymbol{J}_{D}$ in the frequency domain (Section VI).

• •

Provide an algorithm to solve the transceiver design problem (18) by alternating between optimizing the transmit $g(t)$ and the receive filter $h(t)$ (Section VIII).

• •

Explore the Pareto-optimal region obtained by jointly optimizing the transmit and receive filter $g(t)$ and $h(t)$ with respect to a weighted MSE criterion (Section X).

• •

Discuss the ambiguity function resulting from the optimized transceiver design (Section X).

• •

Verify the performance obtained under an optimized design by Monte-Carlo simulation of an asymptotically efficient algorithm $\boldsymbol{\hat{\theta}}(\boldsymbol{y})$ (Section XI).

Our final conclusions are outlined in Section XII. Note that preliminary results have been discussed in [6, 7, 8], where we have focused on either optimizing the receive or the transmit filter of a sub-Nyquist system, and [9] which was centered around a compressive sensing framework.

IV Related Work

The analysis and optimization of wireless transceiver systems enjoys significant attention in the signal processing community [10, 11, 12, 13]. While for example [14, 15] focus on signal optimization for classical radar systems, waveform design of orthogonal frequency division multiplexing (OFDM) signals for radar applications have been considered in [16, 17]. The favorable design of satellite-based positioning signals with high time-delay estimation accuracy is discussed in [18]. Optimization of transmit and receive filters in the spatial domain plays a crucial role in multiple-input multiple-output (MIMO) radar systems [19, 20]. For cognitive radars, high accuracy is achieved by using a Bayesian tracking algorithm and adapting the transmit signal to the state-space knowledge [21]. In [22] the concept of cognitive radar is extended to MIMO radar systems resulting in improved parameter identifiability and higher resolution [23] in comparison to phased arrays. Further, MIMO radar features the effect of virtual receive antennas [24], providing high accuracy at lower hardware costs.

When performing the transceiver system optimization various figures of merit have been considered. In particular mutual information and signal-to-noise ratio (SNR) are popular metrics [12, 14, 15, 25], while a low MSE for the estimated parameters is an alternative choice and leads to the use of estimation theoretic error bounds [26, 27, 22]. Further, the ambiguity function (AF) of a given transmit signal is a classical tool to characterize the waveform quality for applications where joint estimation of time and frequency shifts is desired [28, 29, 30, 19, 31, 32]. Subspace-based techniques for delay-Doppler estimation problems are presented in [33, 34] and lead to significant complexity reduction of the maximum-likelihood (ML) estimator.

An important aspect in the context of this article is the field of compressed sensing (CS) and finite rate of innovation (FRI) techniques, where sparse signal structures [35] and signals with finite degrees of freedom [36] are exploited to facilitate signal processing. CS-related work focuses on methods for reconstructing signals that embody sparsity from a small number of samples [37], while FRI techniques attempt to reconstruct the original signal without bandwidth limitations [38]. In [39, 40, 41] reconstruction of FRI signals in noise is considered. The work of [42] describes CS methods for parameter estimation applications where perfect signal reconstruction is not required. Taking into account sampling hardware complexity, Xampling has been introduced as a sub-Nyquist sampling technique on multiple analog projections of the received signal [43, 44]. With such an approach, low-rate analog-to-digital conversion architectures for radar systems have been proposed which exhibit the same estimation performance as receivers that operate at the Nyquist rate when the SNR is sufficiently high [45, 46]. In [47] an information-theoretic method is presented to optimize the analog projection functions for compressed sub-Nyquist sampling within a Bayesian estimation framework. It is shown that the MSE of the resulting time-delay estimate can be improved for the undersampled system compared with Nyquist sampling methods. In contrast to these works, we consider likelihood-oriented signal processing and transceiver optimization while using a conventional baseband receiver which operates at a sampling rate $f_{s}$ smaller than the Nyquist bandwidth $B$.

V Delay-Doppler Estimation

V-A Specific Receive Model

To illustrate a specific case of the generic transceiver design problem (18), we consider a single-path propagation model,

[TABLE]

with path gain $\gamma\in\mathbb{C}$, time-delay $\tau\in\mathbb{R}$ and Doppler shift $\nu\in\mathbb{R}$. Such a receive model is often employed in radar systems, where the signal is delayed according to the distance of the reflector and frequency shifted (Doppler effect) due to its relative velocity. For solving the filter design problem (18), we ignore $\gamma$ and treat it as known and deterministic, so that the unknown random parameter vector is

[TABLE]

Several reasons motivate this approach. First, accurate estimation of the linear coefficient $\gamma$ is less crucial for most applications where the delay and Doppler are used for ranging and localization, and thus we focus on the optimization of the estimation accuracy for these parameters. This leads to a lower BCRLB for $\tau$ and $\nu$ than for the case where $\gamma$ is treated as unknown [3, p. 48]. The resulting BCRLB still serves as a valid lower bound for the estimation error of $\tau$ and $\nu$. Second, considering $\gamma$ as known simplifies the optimization of the BCRLB and reduces the computational effort when evaluating the corresponding Pareto-optimal region. In practice, $\gamma$ is a nuisance parameter that must still be estimated. We will illustrate in Section XI that the MSE for $\tau$ and $\nu$ is still very close to the smaller BCRLB even when $\gamma$ must be estimated. This indicates that there is little loss when ignoring $\gamma$ during the optimization of the transceive filters $g(t)$ and $h(t)$.

In the following it is assumed that the unknown parameters are Gaussian distributed

[TABLE]

and stochastically independent, such that the covariance matrix of the prior is

[TABLE]

V-B Frequency Domain Representation

Solving the optimization problem (18) requires an analytical characterization of the EFIM (11) under the receive model (19). A frequency-domain representation enables a compact and insightful notation for the received signal model [7, 48] and thus of the FIM entries (13).

For the derivation, we assume a fixed sampling rate $f_{s}$ at the receiver. The transmit filter $g(t)$ is band-limited with two-sided bandwidth $B$. In contrast to the common assumption $B\leq f_{s}$, in our setup we allow $B>f_{s}$. Due to periodicity, the waveform $\breve{x}(t)$ can be represented by its Fourier series

[TABLE]

where $\omega_{0}=\frac{2\pi}{T_{0}}=2\pi f_{0}$. Note that $K=\lceil\frac{2\pi B}{\omega_{0}}\rceil\in\mathbb{N}$ is the total number of harmonics. $G_{k}$ denotes the $k$-th Fourier coefficient of the transmit filter

[TABLE]

We write $\operatorname{\mathcal{K}}$ for the set of Fourier coefficient indices

[TABLE]

with cardinality $|\operatorname{\mathcal{K}}|=K$. Inserting expression (23) into (19) and applying the filtering operation in (3), we obtain

[TABLE]

where

[TABLE]

is the Fourier transform of the receive filter $h(t)$. Evaluating $v(t;\boldsymbol{\theta})$ at sampling instants $nT_{s},n\in\{-\frac{N}{2},\dots,\frac{N}{2}-1\}$ yields

[TABLE]

The set $\operatorname{\mathcal{K}}$ can be partitioned into $N$ disjoint subsets $\mathcal{F}_{k},k\in\{-\frac{N}{2},\dots,\frac{N}{2}-1\}$, where each $\mathcal{F}_{k}$ contains elements of $\operatorname{\mathcal{K}}$ in intervals of $N$ and $k\in\mathcal{F}_{k}$. Denote by $L^{-}_{k}$ the integer such that $k-NL^{-}_{k}$ is the smallest element of $\mathcal{F}_{k}$, and respectively $L^{+}_{k}$ the integer such that $k+NL^{+}_{k}$ is the largest element of the subset $\mathcal{F}_{k}$. The subsets $\mathcal{F}_{k}$ then can be denoted as

[TABLE]

Splitting the sum over $k$ in (28) according to this partition and further using $\omega_{s}=2\pi f_{s}=N\omega_{0}$ yields

[TABLE]

Here $L^{-}_{k}$ and $L^{+}_{k}$ can be interpreted as the number of signal spectrum points in the negative, respectively positive frequency domain that superpose at the frequency point $k\omega_{0}$ due to sampling at a rate $\omega_{s}$. These limits are given by

[TABLE]

Note that if the bandwidth is an odd multiple of the sampling frequency $f_{s}$, i.e,

[TABLE]

the number of overlapping frequencies in the positive and negative frequency domain adjoin and do not depend on the center $k$. In this case $L^{-}_{k}$ and $L^{+}_{k}$ are

[TABLE]

Stacking the entries (V-B) into a vector yields

[TABLE]

Denoting $\circ$ as the element-wise Hadamard product, (35) contains an aliasing matrix $\boldsymbol{A}\in\{0,1\}^{N\times K}$ with entries

[TABLE]

a diagonal time-delay matrix $\boldsymbol{T}(\tau)\in\mathbb{C}^{K\times K}$ with components

[TABLE]

and the transmit filter spectrum vector $\boldsymbol{\tilde{g}}\in\mathbb{C}^{K}$ with elements

[TABLE]

Throughout this work vector and matrix entries are referenced by positive integers, starting at $1$ and therefore $i\in\{1,2,\dots,N\},j\in\{1,2,\dots,K\}$ in equation (38) and $i\in\{1,2,\dots,K\}$ in both preceding equations. The vector $\boldsymbol{\tilde{h}}(\nu)\in\mathbb{C}^{K}$ in (35) denotes the frequency shifted receive filter spectrum with entries

[TABLE]

Further, in (35) we write $\boldsymbol{W}\in\mathbb{C}^{N\times N}$ for the DFT matrix with entries

[TABLE]

where $i,j\in\{1,2,\dots,N\}$ and $\boldsymbol{\Delta}(\nu)\in\mathbb{C}^{N\times N}$ represents the diagonal matrix with elements

[TABLE]

Introducing the Doppler convolution matrix

[TABLE]

a signal representation in the frequency domain

[TABLE]

is obtained from (35), where

[TABLE]

is the attenuated and delay-Doppler shifted received signal spectrum. A frequency domain representation of the noise covariance matrix (6) is obtained by the transformation

[TABLE]

Note that due to the Doppler shift in (41), frequencies $k\omega_{0}+2\pi\nu$ of the receive filter $H(\omega)$ have to be considered for evaluation of the filtered and sampled signal $\boldsymbol{v}(\boldsymbol{\theta})$. The entries of the FIM $\boldsymbol{J}_{F}(\boldsymbol{\theta})$ are directly obtained by plugging (V-B) into (13) and can be found in Appendix -A.

VI Approximation of the Delay-Doppler FIM

These expressions for the FIM make it possible to formulate the transceiver optimization problem (18) for the delay-Doppler model (19). Due to the assumption of a periodic transmit sequence $\breve{x}(t)$, the minimization over the transmit filter $g(t)$ reduces to an optimization with respect to the $K$ entries of $\boldsymbol{\tilde{g}}$. However, using the signal characterization (V-B) does not lead to an analytically tractable optimization problem, as it is a function of the filter spectrum $H(\omega)$ evaluated at frequencies $k\omega_{0}+2\pi\nu$ (see, e.g., (28)). By the definition of the EFIM (11) this leads to an objective function which is the expected value of a function that depends on $H(k\omega+2\pi\nu)$ with continuous argument $\nu$. Therefore in the following section we discuss several steps to approximate the FIM entries and formulate a tractable optimization problem (18). The aim of this section is to simplify the optimization problem so that it can be solved analytically. We then illustrate in Section X that the optimization over the approximated EFIM leads to a significant gain for the exact EFIM. We further show that the region spanned by the approximated EFIM and the exact EFIM exhibit close correspondence. In Section XI, we finally show that for the considered delay-Doppler estimation problem the BCRLB, which by definition involves the exact EFIM, accurately characterizes the performance of the MAP estimator.

VI-A FIM Approximation - Periodic Doppler

To decouple the receive filter variable $\boldsymbol{\tilde{h}}(\nu)$ and the propagation parameter $\nu$, we approximate the Doppler shift to be periodic with $T_{0}$, i.e.,

[TABLE]

where we use a centered version of the modulo operation

[TABLE]

Note that in the interval t\in\big{[}-\frac{T_{0}}{2},\frac{T_{0}}{2}\big{)} the approximation (48) is exact. However, due to the filtering operation $h(t)$ portions of the signal outside that time range are convolved into the sampling interval and thus produce an approximation error.

Under the assumption (48), it is possible to write

[TABLE]

Calculating the Fourier coefficients $d_{z}(\nu)$ yields

[TABLE]

Replacing the Doppler shift in (19) with (48) yields the signal approximation

[TABLE]

with the substitution $z=m-k$. Using the abbreviation

[TABLE]

the signal samples are given by

[TABLE]

Neither the time-delay nor the Doppler shift have a significant impact on the signal bandwidth $B$. The time-delay merely adds a linear phase to the signal spectrum and the Doppler stretch shifts the spectrum by $2\pi\nu$ (see, e.g., (28)). Since $2\pi\nu\ll\omega_{0}$ for practical applications, the receive filter can be assumed to have bandwidth $B$. It follows that $H_{k}=0$ for all $k$ with $|k\omega_{0}|>\pi B$. Hence,

[TABLE]

In vector notation, we obtain

[TABLE]

where the approximated received signal spectrum is

[TABLE]

with the Doppler convolution matrix

[TABLE]

and the receive filter spectrum $\boldsymbol{\tilde{h}}\in\operatorname{\mathbb{C}}^{K}$

[TABLE]

The advantage of the periodic Doppler approximation (48) is that the filter spectrum $\boldsymbol{\tilde{h}}$ does not depend on the Doppler $\nu$. This implies that due to the periodicity of $\bar{v}(t;\boldsymbol{\theta})$ only the harmonics $k\omega_{0}$ of the receive filter $H(\omega)$ are excited and other parts of the filter spectrum do not impact the description of the received signal. Defining the matrix

[TABLE]

the received signal spectrum can be represented by

[TABLE]

VI-B FIM Approximation - Circular Noise Covariance

Computation of the FIM entries (13) requires the inverse of the noise covariance matrix $\boldsymbol{R}_{\boldsymbol{\eta}}^{-1}$. The fact that $\boldsymbol{R}_{\boldsymbol{\eta}}$ is a non-trivial function of the receive filter $h(t)$ hinders tackling the optimization problem (18). We exploit the fact that the noise covariance $\boldsymbol{R}_{\boldsymbol{\eta}}$ is approximately diagonalized by the DFT matrix for large $N$ [48], i.e.,

[TABLE]

where $\boldsymbol{\Omega}_{\boldsymbol{\eta}}\in\operatorname{\mathbb{R}}^{N\times N}$ is a diagonal matrix with entries

[TABLE]

Note that (66) results from Szegös theorem [49] applied to Toeplitz matrices. Defining the receive filter spectrum matrix

[TABLE]

we directly find the useful expression

[TABLE]

VI-C Approximated FIM Entries

With the approximations (48) and (66), the FIM elements (12) are directly obtained by inserting (56) and (66) into (13)

[TABLE]

where $\boldsymbol{\bar{J}}_{F}(\boldsymbol{\theta})$ denotes the FIM that stems from the periodic Doppler and noise covariance approximation. Using the diagonalization property for the Hadamard product of two vectors, equation (70) can be denoted by

[TABLE]

with

[TABLE]

The FIM expression (71) allows us to approximate the objective function (18) by

[TABLE]

with

[TABLE]

Alternatively to (71), it is possible to emphasize the dependence of the FIM on $\boldsymbol{\tilde{h}}$ and reformulate (70) as

[TABLE]

where

[TABLE]

The filter spectrum points

[TABLE]

separated by $\omega_{s}$, experience aliasing since they are combined by the multiplication of $\boldsymbol{\xi}_{i}(\boldsymbol{\theta})\circ\boldsymbol{\tilde{h}}$ with $\boldsymbol{A}$, where the corresponding coefficients are

[TABLE]

Thus, the compact FIM representation

[TABLE]

is obtained. Note that in (78) the summand $\frac{K}{2}+1$ assures a correct indexing of the vector $\boldsymbol{\xi}_{i}(\boldsymbol{\theta})$. In contrast to (VI-C), eq. (79) allows us to characterize the objective function by

[TABLE]

with

[TABLE]

Note that (VI-C) emphasizes the quadratic dependence of the cost function $\mathrm{tr}(\boldsymbol{M}^{\prime}\boldsymbol{\bar{J}}_{D})$ on the transmit spectrum $\boldsymbol{\tilde{g}}$, while (80) is written as a function of the receive filter $\boldsymbol{\tilde{h}}$.

VII Symmetry of the EFIM

The EFIM exhibits symmetry concerning time and frequency reversal of the two filters $\boldsymbol{\tilde{g}}$ and $\boldsymbol{\tilde{h}}$. We will see later that this observation can be used to accelerate the optimization algorithm solving the transceiver design problem (18).

VII-A EFIM with Mirrored Transceive Filter Spectra

Define the mirroring matrix

[TABLE]

which obeys

[TABLE]

such that the frequency mirroring operations $X(-\omega)$ and $H(-\omega)$ can be denoted as $\boldsymbol{\Pi}\boldsymbol{\tilde{g}}$ and $\boldsymbol{\Pi}\boldsymbol{H}\boldsymbol{\Pi}$, respectively. Writing \boldsymbol{\bar{J}}_{F}(\boldsymbol{\theta})\big{|}_{\boldsymbol{\tilde{g}},\boldsymbol{\tilde{h}}} to denote the FIM with transmit $\boldsymbol{\tilde{g}}$ and receive filter spectrum $\boldsymbol{\tilde{h}}$, we find the FIM entries after mirroring the transceive filters along the frequency axis

[TABLE]

It can be seen that (VII-A) contains the mirrored matrix

[TABLE]

With the structure of the time-delay matrix (39), we have

[TABLE]

Note that by definition (39) the matrix $\boldsymbol{T}(\tau)$ is not perfectly symmetric as the frequency interval ranges from $-\frac{K}{2}\omega_{0}$ to $(\frac{K}{2}-1)\omega_{0}$. However, for large $K$ the impact of the single asymmetric frequency point $-\frac{K}{2}\omega_{0}$ is negligible. Using the symmetry of the Doppler coefficients

[TABLE]

we obtain

[TABLE]

Thus,

[TABLE]

showing that applying the mirroring operation to both filters corresponds to a reversal of the unknown parameters, i.e.,

[TABLE]

The time-delay and Doppler shift are assumed to be normally distributed and thus have symmetric distributions, i.e.,

[TABLE]

In general, mirroring a function when averaging with respect to a symmetric distribution does not influence the expected value and thus we obtain the approximate identity

[TABLE]

Note, that using the properties of the Fourier transform, the same result is obtained when mirroring the transmit signal and received signal in the time domain.

VII-B Off-Diagonal Elements of the EFIM with Symmetric Filters

We will show in the following that for symmetric transceiver filters the non-diagonal elements of the FIM vanish. Due to the symmetry and by the properties of the Fourier transform, $\boldsymbol{\tilde{g}},\boldsymbol{\tilde{h}}$ become real-valued. Therefore, we have

[TABLE]

It is straightforward to show that the following properties hold for the complex conjugates

[TABLE]

For the off-diagonal elements of the FIM, the matrices $\boldsymbol{\Phi}_{ij}(\boldsymbol{\theta})$ contain exactly one partial derivative with respect to $\tau$ and one with respect to $\nu$. Hence, for $i\neq j$

[TABLE]

Due to the symmetry of the time-delay distribution $p(\tau)=p(-\tau)$, taking the expected value of (96), the off-diagonal elements of the EFIM $\boldsymbol{\bar{J}}_{D}$ become

[TABLE]

VII-C Compact Cost Functions

With (92), the reversal of the transmit and the receive filter in the frequency domain does not significantly change the elements of the EFIM if the parameters are distributed symmetrically. Therefore, we will simplify the optimization problem (18) by exclusively considering transceive filters $\boldsymbol{\tilde{g}}$ and $\boldsymbol{\tilde{h}}$ which have symmetric impulse responses and spectra. Note, that due to the properties of the Fourier transform and filter symmetry, both filters are real-valued and optimization is performed by adapting half of the spectrum.

We split the transmit signal

[TABLE]

with the one-sided transmit spectrum $\boldsymbol{\tilde{g}}_{\mathrm{r}}\in\operatorname{\mathbb{R}}^{\frac{K}{2}-1}$ and $\boldsymbol{\tilde{g}}_{\mathrm{f},\mathrm{r}}=\boldsymbol{\Pi}\boldsymbol{\tilde{g}}_{\mathrm{r}}$. Note that the first element in (98) has been set to zero, i.e., $G_{-\frac{K}{2}}=0$, to obtain a fully symmetric description. This allows us to rewrite the cost function (VI-C) as

[TABLE]

with the matrices and vectors

[TABLE]

and a reduced filter matrix

[TABLE]

It can be seen that with symmetric filters the cost function (VI-C) can be written in a compact quadratic form (115), requiring only one side of the transmit filter spectrum $\boldsymbol{\tilde{g}}$. A compact version of (80) can be found in a similar way by using

[TABLE]

This allows to rewrite the cost function (80) as

[TABLE]

VIII Transceiver Optimization Algorithm

In the following we focus on solving the transceiver design problem (18) with the approximate EFIM derived in the previous section. Further, we will present an algorithm which solves the maximization task (121). With the frequency domain representation (40) and (63), the optimization problem (18) can be stated as

[TABLE]

Note that in contrast to (18), where the exact EFIM is used, the problem is now formulated based on the EFIM approximation $\boldsymbol{\bar{J}}_{D}$ that has been introduced in Section VI.

We first focus on the problem of optimizing the transmit signal $\boldsymbol{\tilde{g}}$ for a given receive filter $\boldsymbol{\tilde{h}}$

[TABLE]

where the cost function is expressed by (VI-C). It is well known that the solution to such a problem is the Eigenvector $\boldsymbol{\gamma}_{1}$ of the matrix $\boldsymbol{\Phi}$ corresponding to the largest Eigenvalue [50].

On the other hand, for a constant transmit signal $\boldsymbol{\tilde{g}}$, the optimization problem reads as

[TABLE]

Having a closer look at the filter representation of the objective function (80), it is straightforward to prove that two different vectors $\boldsymbol{\tilde{h}}_{k_{1}}$ and $\boldsymbol{\tilde{h}}_{k_{2}}$, $k_{1}\neq k_{2}$ consist of two disjoint sets of filter coefficients $H_{k}$. Therefore problem (123) is solved by independently maximizing each term of the sum (VII-C), which results in choosing $\boldsymbol{\tilde{h}}_{k}=\boldsymbol{\delta}_{k,1}$ where $\boldsymbol{\delta}_{k,1}$ is the Eigenvector of $\boldsymbol{\Delta}_{k}$ corresponding to the largest Eigenvalue.

With the two separate optimization tasks (122) and (123) for the transmit and the receive filter, we suggest the following iterative algorithm to find the solution to the design problem in (121). The algorithm is initialized with filters $\boldsymbol{\tilde{g}}_{\mathrm{init}}$ and $\boldsymbol{\tilde{h}}_{\mathrm{init}}$ that can either be chosen randomly or set to initial prototype versions. In an alternating way, in the $i$-th iteration the analog transmit filter $\boldsymbol{\tilde{g}}$ is optimized by solving (122) for a fixed receive filter $\boldsymbol{\tilde{h}}^{i-1}$. The solution $\boldsymbol{\tilde{g}}^{i}$ is then used to improve the solution for $\boldsymbol{\tilde{h}}$ with (123). This approach ensures that the value of the objective function increases with each iteration. If the relative performance gain is less than a threshold value $\epsilon$, the algorithm terminates and outputs the solutions $\boldsymbol{\tilde{g}}^{\star}$ and $\boldsymbol{\tilde{h}}^{\star}$. Note that such an alternating approach does not necessarily find the globally optimum solution. However, simulations show (see Sections X and XI) that the procedure outputs a transceiver design which provides significant performance gains. For all tested scenarios, the suggested algorithm converges quickly and returns a filter pair with favorable design after a few steps.

IX Reference Systems

To provide a comparison to commonly used transceive filters, we introduce two reference designs which will serve as a benchmark for the optimized solutions. We use rectangular pulses as employed in Global Navigation Satellite Systems (GNSS) [51] such as the American Global Positioning System (GPS), and linear frequency modulation (LFM) pulses that find application in long-range, high-resolution radar systems [52, p. 88]. For further details about LFM pulses, the reader is referred to [53].

Both reference signals $\breve{x}(t)$ can be described by the periodic repetition of a transmit pulse $g(t)$ in intervals of $T_{0}$. Having available the spectrum $G(\omega)$ of the reference pulse, the transmit filter Fourier coefficients are simply obtained by the evaluation of the discrete spectrum points $G(k\omega_{0})$ and normalizing to $\boldsymbol{\tilde{g}}^{\mathrm{H}}_{\mathrm{ref}}\boldsymbol{\tilde{g}}_{\mathrm{ref}}=P_{T}$. The analog reference receive filter $\boldsymbol{\tilde{h}}_{\mathrm{ref}}$ is an ideal low-pass filter with bandwidth $B_{\mathrm{ref}}=f_{s}$ which complies with the sampling theorem. To visualize the possible performance gains, in the following, we consider the relative performance measures

[TABLE]

which characterize the information gain in the data influenced part of the BCRLB (9) and therefore characterize the factor by which the transmit power of the optimized system could be reduced while maintaining the performance level of the reference system. Note that (124) is computed using the exact EFIM (11) with entries (153) – (155). The reference systems are configured such that they use the same transmit power $P_{T}$ as the optimized system.

IX-A Rectangular Phase Code Modulation

The rectangular prototype pulse has duration $T=2T_{s}$ and its Fourier transform can be denoted

[TABLE]

Under the assumption that the pulse is ideally band-limited to a two-sided bandwidth $B=f_{s}$ at the transmitter, we obtain

[TABLE]

The rectangular phase code (RPC) is then generated by convolving the prototype rectangular pulse with a code sequence

[TABLE]

where

[TABLE]

for some known binary sequence $\beta_{m}\in\{-1,1\}$.

IX-B Linear Frequency Modulation (LFM)

The LFM reference signal consists of an envelope $a_{\mathrm{LFM}}(t)$ and a sinusoid with linearly increasing frequency $w_{\mathrm{LFM}}(t)$. For the case of a rectangular envelope, we have

[TABLE]

with pulse duration $T$ and frequency slope $\mu$. In the literature, (134) is often referred to as chirp pulse. The frequency of the sinusoid is increasing linearly in time

[TABLE]

The LFM transmit pulse and its spectrum are [52, p. 138]

[TABLE]

with the complex Fresnel integral $Z(u)=\int_{0}^{u}\mathrm{e}^{\mathrm{j}\pi\frac{a^{2}}{2}}\mathrm{d}a$. The two-sided bandwidth of this pulse is approximately $B_{\mathrm{LFM}}=\frac{\mu T}{2\pi}$. For the reference system, we use a pulse length of $T=T_{0}$. The bandwidth $B_{\mathrm{LFM}}=\frac{2f_{s}}{3}$ uses the full frequency band without violating the sampling theorem and yields $\mu=\frac{4\pi N}{3T_{0}^{2}}$.

X Optimization Results

X-A Performance - Pareto-optimal Transceiver Region

In general there exists a trade-off between the estimation accuracy of the time-delay and the Doppler shift. As such, for the presented transceiver optimization framework we examine the Pareto-optimal set $\mathcal{P}=\{(g(t),h(t))\}$ of transmit and receive filters, for which the estimation of one parameter cannot be improved through the filter design without reducing the accuracy of the other parameter. We define the set $\mathcal{G}$ of relative performance gains (124) which are obtained by filters on the boundary of the Pareto-optimal region $\mathcal{P}$

[TABLE]

The set $\mathcal{G}$ can be approximated by solving the transceiver optimization problem (121) over all positive definite weighting matrices $\boldsymbol{M}^{\prime}$ and computing the corresponding EFIM values (11). In Section VII-B it has been shown that for the examined system the non-diagonal elements of the EFIM $\boldsymbol{\bar{J}}_{D}$ vanish such that it is sufficient to consider the diagonal weighting matrices

[TABLE]

with $\alpha\in[0,1]$.

X-A1 Variable Bandwidth with Fixed Sampling Rate

Fig. 2 shows the Pareto boundary for different bandwidths $B=(2L+1)f_{s}$ and signal periods $T_{0}$, where $f_{s}=25\mathrm{MHz}$ has been used. Here, the rectangular transmit waveform (131) with $T_{0}=2\mu\mathrm{s}$ is used as the reference in (124) and the variances of the random unknown parameters (22) are set to $\sigma_{\tau}=1\mathrm{ns},\sigma_{\nu}=5\mathrm{kHz}$. Fig. 2 has been created by solving (18) using the approximated FIM (70) for symmetric signals under all weightings $\alpha$ and computing the corresponding exact values $\chi_{\tau/\nu}$ (124). Note that the system $L=0,T_{0}=4\mu\mathrm{s}$ operates at half the transmit power $\frac{P_{T}}{2}$ such that the total transmitted energy $P_{T}T_{0}$ is the same as in the reference system. It can be observed that for the estimation of the time-delay $\tau$, increasing the bandwidth $B$ leads to a substantial performance gain, while for the estimation of the Doppler, an extension of the time-spread $T_{0}$ improves the achievable accuracy. The fact that the Pareto region appears to be larger for the Doppler parameter when the bandwidth is lower can be explained by the inaccuracies that result from the approximation of the EFIM used in (121). We see that for essentially any choice of the weighting matrix, the optimized design provides significant gains in estimation performance compared with the prototypes. We also note that the nearly rectangular Pareto regions indicate that, with optimized filter designs, the trade-off between delay and Doppler estimation performance is relatively small.

X-A2 Fixed Bandwidth with Variable Sampling Rate

We have seen that optimized filters have the potential to increase the estimation performance of delay-Doppler estimation methods. We now investigate the estimation performance for a fixed transmit bandwidth $B=25\mathrm{MHz}$, a signal period $T_{0}=2.4\mu\mathrm{s}$ and different sampling frequencies $f_{s}=\frac{B}{2L+1}$.

Fig. 3 shows the Pareto regions of the optimized waveforms with respect to a rectangular signal, all with the same bandwidth of $B=25\mathrm{MHz}$. Hence, the number of Fourier coefficients $K=BT_{0}$ is the same for all systems. Note that the sampling rate for the reference system is held constant, while the sampling rate and the number of samples $N$ of the optimized system decrease with increasing $L$. This indicates that although lower sampling rates are used, the optimized waveform design can still provide high estimation accuracy. In particular, this means that, if the bandwidth $B$ is fixed and the transceive filters are chosen according to the proposed design rule, reducing the number of samples inside the constant time interval $T_{0}$ does not result in an accuracy loss for the estimation of $\tau$. There are two main reasons that explain this phenomenon. First, the time delay estimation accuracy benefits from transmit spectra whose power is concentrated at high frequencies (see, e.g., Fig. 5a). By violating the design rule $B\leq f_{s}$, the transmit spectrum can be chosen to feature bandwidths with high frequencies, although the sampling bandwidth is decreased. Second, the power loss induced by the reduced number of samples is compensated by a bandpass structure of the receive filter (see Fig. 6d), which decreases the power of the noise accordingly. However, the Doppler estimation suffers slightly from a reduced number of sampling points.

X-A3 Comparison of Approximate and Exact EFIM

In the previous paragraphs, the transmit and receive filter optimization was based on the approximated EFIM. In this section, we investigate the accuracy of the noise covariance and Doppler approximation. We show that evaluating the approximated and exact EFIM for the optimized systems leads to similar Pareto regions.

Figure 4 shows the Pareto-optimal region of the approximated and exact Fisher information matrix, for the case of a rectangular reference signal (131). Fig. 4 was created by solving the transceiver optimization problem (18) using the approximated EFIM (VI-C) under all weightings $\alpha$, as illustrated in Section VIII. The optimized filters were then used to compute the corresponding exact performance gains $\chi_{\tau/\nu}$ (124) using (153)-(155). For comparison, approximated values $\bar{\chi}_{\tau/\nu}$ have been calculated using (70). The system parameters were chosen equal to the earlier simulations, i.e., $f_{s}=25\mathrm{MHz}$ and $B=(2L+1)f_{s}$. The parameter variances are $\sigma_{\tau}=1\mathrm{ns}$, $\sigma_{\nu}=5\mathrm{kHz}$. It can be observed that the approximated bound lies close to the exact bound, with a gap of maximum $2\mathrm{dB}$ and exhibits a similar shape. Note that the optimization in Section VIII is performed by using the approximated information measure (VI-C). Therefore higher gains might be achievable when using a more accurate representation of the exact EFIM.

X-B Design - Optimized Transmit and Receive Filters

To illustrate transceiver designs optimized for the estimation of $\tau$ and $\nu$, we visualize time and frequency domain representations of the analog receive and transmit filters obtained with the presented optimization algorithm. For the same configuration as in Fig. 2 and a signal period of $T_{0}=2\mu\mathrm{s}$ we first investigate the case where the sampling theorem is satisfied, i.e., $B=f_{s}$. Note that under these conditions the optimum receive filter is an ideal rectangular low-pass filter with two-sided bandwidth $B_{\mathrm{LP}}=f_{s}$.

Fig. 5a shows the spectrum of the optimized transmit signal in the frequency-domain when full weight is given to the optimization of the delay estimation ($\alpha=1$).

In contrast Fig. 5b depicts the result of the transmit filter optimization in the time-domain when the inference of the Doppler is emphasized ($\alpha=0$). It can be observed that when focusing on the estimation of the delay ($\alpha=1$), the transmit signal power is concentrated at high frequencies (see Fig. 5a). When enforcing estimation of the Doppler, the signal power is mainly located at the boundaries of the observation interval (see Fig. 5b). This is because for delay estimation, signal components with high frequencies are more beneficial for the EFIM, since the phase shift, which is introduced by the propagation delay, increases linearly in frequency. On the other hand, for Doppler estimation a signal with wide time-spread is beneficial as the phase shift that is introduced by the Doppler effect increases linearly in time. Using a weighting

[TABLE]

which simultaneously emphasizes both the estimation of the delay and Doppler shift, we obtain the transmit filter depicted in Figs. 5c and 5d. Note that $\alpha^{*}$ leads to the point on $\mathcal{G}$ at which the sum of both accuracy gains (124) is maximum, i.e., the most upper right point in the Pareto region $\mathcal{G}$. It can be observed that in this case, the signal power is concentrated at the boundaries of the time and frequency domain. Due to the properties of the Fourier transform, the signal cannot be limited in time and frequency simultaneously, which leads to the observed trade-off between delay and Doppler estimation. Furthermore, it is noticeable that the optimized transmit signal exhibits strong oscillations in the time and frequency domain (Figs. 5c and 5d).

Allowing for a higher bandwidth $L=1$, i.e., $B=3f_{s}$ and therefore for received signals which violate the sampling theorem, Figs. 6a and 6b show the optimized waveform for the max-sum weighting $\alpha^{\star}$ (140). The dashed lines visualize the interval $[-\frac{\omega_{s}}{2},\frac{\omega_{s}}{2})$. It can be observed that the optimized transmit filter has a total bandwidth of $B>f_{s}$ and therefore violates the sampling theorem. However, as this signal shows a bandpass character, the real bandwidth of the signal, i.e., the cumulative portion of the frequencies where the signal exhibits non-zero spectral power density is merely $\omega_{s}$. Although the total bandwidth of the signal exceeds $\omega_{s}$, due to its bandpass character the transmit signal does not experience significant aliasing during the sampling process at the receiver. Comparing the sampled signal spectrum in Fig. 6c with the transmit spectrum in Fig. 6a, it can be seen that signal replicas are not distorted. The advantage of the band-pass character of the transmitter is a high sensitivity at the receiver with respect to the propagation delay. Fig. 6d displays the optimized receive filter spectrum, which is a bandpass. An interesting observation is that none of the optimized signals uses the full available bandwidth $B$. Note that due to the relatively small prior variances, in general, the full bandwidth and time interval could be used. For higher parameter uncertainties, safety margins at the interval boundaries would occur to prevent a potential power loss from large time or frequency shifts.

X-C Complexity - Ambiguity Function

The ambiguity function, defined as the output of a matched filter without Doppler shift compensation

[TABLE]

is a well-known indicator for analyzing the suitability of transmit waveforms for delay-Doppler estimation[52, ch. 4]. It can be shown that (141) is directly related to the functional distance of the transmit waveform and its time and frequency-shifted versions. It therefore provides insight about the self-similarity of the signal and its sensitivity with respect to shifts in time and frequency. By discussing the curvature of (141) at $\psi(0,0)$ and its side-lobe levels, it is possible to establish a measure of the adequateness of a waveform with respect to delay and Doppler estimation [54].

X-C1 MAP Estimation

In a realistic scenario the attenuation $\gamma$ in the received signal model (19) is deterministic and unknown, so the joint maximum-a-posteriori/maximum-likelihood (JMAP-ML) approach [55]

[TABLE]

has to be used to estimate the unknown parameters. For the case of Gaussian noise, the likelihood function becomes

[TABLE]

Writing

[TABLE]

we find that the solution of (142) needs to satisfy

[TABLE]

Hence, we obtain the estimator for the path gain

[TABLE]

as a function of the parameters $\boldsymbol{\theta}$. Substituting (146) into (142), we obtain the estimator

[TABLE]

as the maximizer of the MAP function

[TABLE]

X-C2 MAP Ambiguity Function

On the basis of (148), we define a function

[TABLE]

closely related to the ambiguity function (141), which we refer to as the MAP ambiguity function (MAP-AF). This measure is derived directly from the likelihood function by an estimation theoretic argument and has the advantage of incorporating the prior knowledge about the unknown parameters $\boldsymbol{\theta}$ into the ambiguity function. Further, in contrast to the standard ambiguity function (141), the uncertainty about the attenuation $\gamma$ as well as the effect of colored noise with covariance $\boldsymbol{R}_{\boldsymbol{\eta}}$ is considered by the suggested ambiguity characterization (149).

Figs. 7a and 7b show the normalized MAP-AF for the optimized transmit and receive filter and for the reference rectangular pulse (131) in a signal-dominant environment, i.e., $N_{0}\rightarrow 0$ to highlight the effect of the transmit and receive filters. It can be seen that the optimized system exhibits a MAP-AF maximum that is significantly narrower than that of the LFM signal which allows for a higher time-delay and Doppler estimation accuracy. Note that the optimized system shows local maxima, which can cause ambiguities for the solution of the estimator (X-C1). However, for a priori knowledge with small uncertainties, i.e., $\sigma_{\nu}\ll f_{0}$, $\sigma_{\tau}\ll T_{s}$, these local extrema do not influence the estimation as the probability that a MAP-AF function side-lobe exceeds the main lobe for a particular noise realization vanishes.

XI Simulation Results

With the estimator (X-C1) it is possible to use Monte-Carlo simulations to verify the achievable performance gain with an optimized transceiver design. To this end, we use $f_{s}=25\mathrm{MHz}$ and $T_{0}=2\mu\mathrm{s}$. The uncertainty of the time and Doppler shift are $\sigma_{\tau}=1\mathrm{ns}$ and $\sigma_{\nu}=5\mathrm{kHz}$ and $L=1$. The optimized transmit and receive filters are obtained by using the max-sum weighting $\alpha^{\star}$ for the matrix $\boldsymbol{M^{\prime}}$. Figs. 8 and 9 show the empirical normalized mean squared error (NMSE)

[TABLE]

of the estimator (X-C1) for the optimized and rectangular reference systems compared with the analytical BCRLB (9) on a double logarithmic scale, where $\mathrm{MSE}_{\hat{\tau}/\hat{\nu}}$ represents the diagonal elements of (8), empirically evaluated based on the estimation in (142).

Note that we use the per-Hz signal-to-noise ratio (pSNR) $\frac{P_{T}}{N_{0}}$, which is independent of the receive filter, to establish a fair comparison between different transceivers. For low pSNR, the BCRLB saturates since the receive data is so noisy that it does not contribute to the performance of the estimator (X-C1). The estimator then exclusively operates on the basis of the prior information and simply produces the mean $\boldsymbol{\hat{\theta}}=\boldsymbol{\mu}_{\boldsymbol{\theta}}$. In the high pSNR regime, it can be observed that the RMSE of the estimator (X-C1) converges to the BCRLB. In this region, the presented transmit and receive filter optimization is beneficial, as the estimator favors the received data in the estimation process. Furthermore, it is observed that for moderate to high pSNR values the performance gain is roughly $20\mathrm{dB}$ for the estimation of the time-delay and $4\mathrm{dB}$ for the Doppler estimation, which corresponds to the findings from the Pareto region in Fig. 2. The pSNR value at which the receive data starts to contribute to the reduction of the estimation uncertainty is approximately $\mathrm{pSNR}_{\tau}=65\mathrm{dB}-\mathrm{Hz}$ for the time-delay and $\mathrm{pSNR}_{\nu}=85\mathrm{dB}-\mathrm{Hz}$ for the Doppler shift. The difference between these two pSNR levels can be explained by the different sensitivity of the transceiver with respect to the two unknown parameters. Note that when using a larger transmit signal duration $T_{0}$, the estimation accuracy of the time-delay increases linearly. In contrast, the Doppler shift performance grows quadratically with $T_{0}$, since the derivative of the receive signal with respect to the Doppler shift scales linearly with the observation interval. As a consequence, the $\mathrm{pSNR}_{\tau}$ and $\mathrm{pSNR}_{\nu}$ threshold level can be aligned by using a larger sampling period $T_{0}$. During our simulations, we observed that the system optimization has no influence on the estimation performance of the linear path gain $\gamma$.

XII Conclusion

We have presented a framework for the joint optimization of the analog transmit and receive filters for estimating the time and Doppler shift with sub-Nyquist sampling. The initial weighted MSE minimization problem has been reformulated as a maximization problem based on the expected Fisher information matrix. Exploiting various approximations, we have shown that it is possible to find transceiver filters that perform significantly better than conventional designs using a simple iterative optimization algorithm. We have demonstrated that the resulting waveforms exhibit bandwidths that violate the sampling theorem ($B>f_{s}$), while the obtained filters have a bandpass structure capturing the entire signal power of the transmit waveform. The Pareto-optimal region exhibits a small trade-off between delay and Doppler estimation and is substantially increased by jointly optimizing the transmit and receive filters. Further, we have modified the standard ambiguity function by taking into account noise covariance and prior information on the delay and Doppler. With the resulting MAP ambiguity function, we have illustrated the gain of the proposed optimized system design relative to conventional approaches. To evaluate the performance of the adapted sub-Nyquist transceiver, Monte-Carlo simulations were performed. The simulation results show significant accuracy gains, such that we conclude that transceiver optimization is an important aspect when designing parameter estimation systems operating under sub-Nyquist conditions.

-A Fisher Information Matrix for Delay-Doppler Estimation

The FIM elements (13) are obtained by computing the derivatives $\frac{\partial}{\partial[\boldsymbol{\theta}]_{i}}\boldsymbol{v}(\boldsymbol{\theta})$. Using (V-B), we obtain

[TABLE]

with the partial derivatives $\boldsymbol{\partial T}(\tau)=\frac{\partial}{\partial\tau}\boldsymbol{T}(\tau)$, $\boldsymbol{\partial\tilde{\Delta}}(\nu)=\frac{\partial}{\partial\nu}\boldsymbol{\tilde{\Delta}}(\nu)$ and $\boldsymbol{\partial\tilde{h}}(\nu)=\frac{\partial}{\partial\nu}\boldsymbol{\tilde{h}}(\nu)$. Inserting (151) and (152) in (13) yields the FIM entries in (153)–(155).

Acknowledgment

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper. They are also grateful to Dr.-Ing. Amine Mezghani for many helpful discussions.