Parameter estimation using the sliding-correlator’s output for wideband propagation channels

Abstract

In this contribution, a high-resolution parameter estimation algorithm is derived based on the Space-Alternating Generalized Expectation-maximization (SAGE) principle for extracting multipath parameters from the output of sliding correlator (SC). The SC allows calculating channel impulse responses with a sampling rate less than that required by Nyquist criterion, and hence is widely used in real-time wideband (e.g., >500 MHz) channel sounding for the fifth generation wireless communication scenarios. However, since the sounding signal needs to be sent repetitively, the SC-based solution is unacceptable for time-variant channel measurements. The algorithm proposed here estimates multipath parameters by using a parametric model of both low- and high-frequency components of the SC’s output. The latter was considered as distortions and discarded in the conventional SC-based channel sounding. The new algorithm allows estimating path parameters with less repetitions of transmitting the sounding signal and still exhibits higher estimation accuracy than the conventional method. Simulations are conducted and illustrate the root mean square estimation errors and the resolution capability of the proposed algorithm with respect to the bandwidth and the length of the SC’s output. These studies pave the way for measuring time-variant wideband propagation channels using SC-based solutions.

1 Introduction

Measurement-based channel models are important for verifying the performance of wireless communication systems in realistic propagation scenarios [1, 2]. Geometry-based stochastic channel models, such as the WINNER spatial channel models [3], IMT-Advanced models [4], and COST2100 multiple-input multiple-output (MIMO) models [5], have been proposed in various standards and widely used to generate single- and multi-link channel realizations at the carrier frequency up to 6 GHz with a bandwidth up to 100 MHz.

Recently, researches on the fifth generation (5G) wireless communications have been paid a lot of attention. The European 7th framework project “Mobile and wireless communications Enablers for the Twenty-twenty Information Society (METIS)” announced a white paper which describes the typical applications and propagation environments considered in 5G [6]. According to the definition by the METIS project, the candidate frequency bands for 5G applications range from 0.45 to 85 GHz, and the bandwidth is from 0.5 up to 2 GHz [6]. At present, the shortage of measurement-based channel models for these frequencies, particularly in the millimeter (mm)-wave bands hinders both the progress of 5G standardization and the designing of 5G-based communication systems and networks. Characterization of mm-wave channel with bandwidth beyond 0.5 GHz for various types of applications and environments began to attract much research attentions recently.

Data acquisition for wideband channels is usually performed by using the equipment such as oscilloscope, spectrum analyzer, and vector network analyzer. The latter two kinds of equipment usually do not have the capability of recording the complex time-domain signals, and thus not suitable for investigating the wideband channel characteristics extracted from multipath parameters. For sampling the mm-wave signals, the oscilloscope devices are required to have sampling rate up to 100 GHz, which is not easy to achieve. Furthermore, due to the small storage in the oscilloscope devices, measurement of wideband channel becomes very time-consuming. A solution tackling these problems is to down-sample the received wideband signals and store the data at a low speed which allows transferring data in real-time from local memory to external disk. Then by using a so-called sliding correlation (SC) technique, a time-dilated approximate channel impulse response (CIR) can be calculated by low-pass filtering (LPF) the received data if the sounding signals can be sent repetitively. The LPF in the receiver is applied to remove the distortion components which have higher frequencies [7]. It has been shown in [8] that pre-filtering techniques can also be applied in the transmitter side to achieve the same objective. Due to the benefits of low complexity in the receiver design and acceptable costs, the SC-based data acquisition has been widely adopted [9–12].

However, the SC-based data acquisition has two problems. First, the higher-frequency components in the SC’s output considered as distortions still carry information of channel parameters and, thus, should be exploited to improve the accuracy of parameter estimation. A drawback resulted when higher-frequency components are considered is that the time-dilated approximation of CIR is unavailable, and as a consequence, conventional peak-searching estimation methods adopted in the SC-based channel estimation are inapplicable. Second, the time-dilated CIR generated by the conventional SC requires the sounding signal being sent repetitively. The number of the repetitions, also called as sliding factor, is usually in the 10³ order of magnitude or even higher [7]. In the cases where channels are time-variant, the CIR may not be calculated within the channel coherence time. As a consequence, the mobile to mobile (M2M) channel measurements cannot be conducted by using the SC-based solution. Recently, a Space-Alternating Generalized Expectation-maximization (SAGE) estimation approach was introduced in [13] which is derived based on a parametric model characterizing the SC’s output, allowing the estimation of multipath parameters by using higher-frequency components. However, this solution still relies on the SC’s output obtained by sending the sounding signals many times. No thorough investigation has been carried out so far for the feasibility of accurate parameter estimation based on the SC’s outputs without sending the sounding signals repetitively.

In this contribution, the SAGE algorithm originally derived in [13] based on a parametric model for both low- and high-frequency components of SC’s output is elaborated. Its performances in estimating multipath parameters are investigated extensively by using simulations. It shows that without discarding the higher-frequency components of SC’s output, the estimation accuracy, particularly for delay parameters, can be improved substantially. In addition, another benefit of this novel estimation algorithm not found previously is discovered; that is, the estimation of path parameters, including Doppler frequency, can be performed by using only a fraction of the SC’s output. Hence, the overall observation span can be kept less than channel coherence time in time-variant cases, and characterizing time-variant channels through SC-based measurements, which cannot be performed before, becomes feasible. Simulations are carried out to compare the performance of the proposed algorithm with the conventional method, and investigate the impact of selecting different bandwidth of LPF and the length of the SC’s output on the RMSEEs, resolution capability of the algorithm.

The rest of the paper is organized as follows. Section 2 describes the parametric signal model. In Section 3, a SAGE algorithm is presented. Section 4 describes the simulation results for the performance of the proposed algorithm. Finally, conclusive remarks are addressed in Section 5. To improve the understandings of the mathematical notations adopted in this contribution, Table 1 lists all the symbols introduced and corresponding explanations.

Table 1 Explanation of adopted symbols

2 Signal model

As elaborated in [9, 10] and [7], the SC performs a specific cross-correlation operation, e.g., between a pseudo-noise (PN) random sequence u(t) with chip rate f _c and another sequence u ^′(t) with chip rate f c′. According to [7], both sequences contain exactly the same chips, and the chip rates are related as $f_{c}'=\frac {\gamma -1}{\gamma }f_{c}$, where γ is called sliding factor. By sample-wise multiplying these two sequences in the time domain for multiple cycles which start with linearly increasing time-offsets and summing the products over individual cycles of u ^′(t), a time-dilated approximate a _u (τ/γ) of the autocorrelation function a _u (τ)=E[u(t)u ^∗(t−τ)] can be calculated by low-pass-filtering the SC’s output with bandwidth B=[−f _c /γ,f _c /γ].

In the channel sounding cases, the received signal is the convolution of the transmitted sequence u(t) with the CIR h(τ), the output of the SC after the LPF with bandwidth B, is the time-dilated approximate $\hat {h}\left (\tau /\gamma \right)$ of the CIR. Here, $\hat {h}\left (\tau /\gamma \right)$ is a time-dilated version of $\hat {h}\left (\tau \right) = h\left (\tau \right)\ast a_{u}\left (\tau \right)$ with ∗ denoting the convolution operation. It is well-accepted that only $\hat {h}(\frac {\tau }{\gamma })$ obtained with the LPF bandwidth B can be used to estimate the channel parameters [7, 14]. Whether the components obtained with larger B, e.g., B _n =[−n f _c /γ,n f _c /γ], n>1, are applicable for estimating the characteristics of h(τ) is necessary to investigate. It is worth mentioning that the time-dilated CIRs with bandwidth of B _n can be obtained by averaging the temporal output of a SC received within the time of L/(f c′n).

Let us consider the case where a time-variant channel consists of M specular paths which are dispersive in the delay and Doppler frequency domains. The baseband representation r(f) of the received signal expressed in frequency domain can be written as

$$ r\left(\,f\right) =\sum\limits_{\ell=1}^{M} \alpha_{\ell} \exp\left\{- j2\pi f\tau_{\ell}\right\} u\left(\,f-\nu_{\ell}\right)+n\left(\,f\right), $$

((1))

where α _ℓ , τ _ℓ , and ν _ℓ denote the complex attenuation, delay, and Doppler frequency of the ℓth propagation path, respectively, and u(f) is the frequency-domain representation of the transmitted maximum-(m-)length PN sequence of L chips. When a rectangular pulse shape is applied, u(f) can be written as [7]

$${} {\fontsize{8.9pt}{9.6pt}\selectfont{\begin{aligned} u\left(\,f\right) = \frac{V_{0}}{L}\sum\limits_{k\in \mathcal{Z}} \text{sinc} \left(\frac{k}{L}\right)\delta\left(f-\frac{f_{c} k}{L}\right)e^{j\frac{k\pi}{L}}\sum\limits_{i=1}^{L}(2 a_{i}-1)e^{-j\frac{2k\pi}{L}i}, \end{aligned}}} $$

((2))

where $\phantom {\dot {i}\!}\mathcal {Z}$ represents the set of integers, V ₀ is the chip magnitude, and a _i ∈[−1,1], i=i,…,L are specified by the m-sequence. The noise component n(f) in (1) is a standard white complex Gaussian random variable:

$$ n\left(\,f\right)= w\left(\,f\right)\sum\limits_{k\in \mathcal{Z}}\delta\left(f-\frac{f_{c} k}{L}\right),\;\;w\left(\,f\right)\sim\mathcal{CN}(0,N_{0}). $$

((3))

The expression u ^′(f) for the lower chip-rate sequence u ^′(t) is similar to u(f) in (2) but with f _c substituted by f c′. The SC’s output y(f)=r(f)∗u ^′(f) is the result of convolution operation in the frequency domain, i.e.,

$$ y\left(\,f\right)=s\left(\,f\right)+n'\left(\,f\right), $$

((4))

where the signal component s(f) can be calculated as

$$ s\left(\,f\right)=\sum\limits_{\ell=1}^{M} \alpha_{\ell} p\left(\,f;\tau_{\ell}, \nu_{\ell}\right) $$

((5))

with

$$p\left(\,f;\tau_{\ell}, \nu_{\ell}\right)=\exp\{- j2\pi f\tau_{\ell}\} u\left(\,f-\nu_{\ell}\right)\ast u'\left(\,f\right) $$

((6))

calculated by invoking the equality δ(f−f ₁)∗δ(f−f ₂)=δ(f−(f ₁+f ₂)) as

$${} \begin{aligned} &p\left(\,f;\tau_{\ell}, \nu_{\ell}\right) =\left(\frac{V_{0}}{L}\right)^{2}\sum\limits_{k,k'\in\mathcal{Z}}e^{-j2\pi \left(\frac{f_{c}k}{L}+\nu_{\ell}\right)\tau_{\ell}-j\frac{\pi}{L}(k+k')}\\&\qquad\qquad\qquad\delta\left(f-\nu_{\ell} - \frac{f_{c}k+f_{c}'k'}{L}\right)\\ &\qquad\qquad\text{sinc}\left(\frac{k}{L}\right)\text{sinc}\left(\frac{k'}{L}\right) \sum\limits_{i=1}^{L}\sum\limits_{i'=1}^{L}\\&\qquad\qquad\qquad\left[(2a_{i}-1)(2a_{i}'-1)e^{-j\frac{\pi}{L}(ki+k'i')}\right] \end{aligned} $$

((7))

and

$$ \begin{aligned} &n'\left(\,f\right)=n\left(\,f\right)\ast u'\left(\,f\right)\\ &=\frac{V_{0}}{L}\sum\limits_{k\in\mathcal{Z}}w\left(\,f\right)\delta\left(\,f-\frac{f_{c}k+f_{c}'k'}{L}\right)\text{sinc}\left(\frac{k}{L}\right)e^{-j\frac{\pi}{L}k'} \\&\quad\sum\limits_{i=1}^{L}(2a_{i}-1)e^{-j\frac{2\pi}{L}(k'i)}. \end{aligned} $$

((8))

The parameters Θ=[α _ℓ ,ν _ℓ ,τ _ℓ ;ℓ=1,…,M] in (4) are unknown and need to be estimated. For simplicity, let us assume that the data y=[y(f);f∈(f ₁,…,f _N )] with N being the total number of frequency bins, obtained within the duration long enough for generating one observation of $\hat {h}(\frac {\tau }{\gamma })$, is available. Estimation of Θ needs to be carried out given y. Notice that this assumption is realistic in the case where the channel coherence time is so short that the SC cannot generate multiple consecutive CIRs for a stationary channel. However, the parameter estimation algorithm derived in the later part of the paper can be easily extended to the case where multiple CIRs are available.

3 Parameter estimation

The maximum likelihood estimate (MLE) of Θ can be derived based on the signal models (4) to (8). However, obtaining the MLE of Θ requires solving a 4M-dimensional optimization problem. The computational complexity involved prohibits any practical implementation. In the following, we present the SAGE algorithm which can iteratively update the subsets of Θ and output the approximate of the MLE of Θ when the estimation process converges [15, 16].

Figure 1 depicts the diagram of the SAGE algorithm derived in the case considered here. To execute the SAGE algorithm, an initial estimate $\hat {\boldsymbol {\Theta }}^{[0]}$ of the unknown parameters Θ is necessary, which can be obtained by using, e.g., the Bartlett beamforming method [17], or parametric approaches based on successive interference cancelation, such as the non-coherent MLE proposed in [18]. The overall parameter estimates are split into multiple subsets. In each iteration of SAGE algorithm, the parameter estimates in a selected subset are updated under the conditions that the observations of received signals are available and that the unknown parameters are assumed to be identical with their estimates that have been calculated from the previous iterations. The SAGE algorithm guarantees that the likelihood of the overall parameter estimates monotonically increases and becomes stabilized after a certain number of iterations, i.e., the so-called “convergence” of the algorithm is achievable. Empirically, when the increment of likelihood as the iteration continues becomes insignificant, or the changes of the parameter estimates compared from a previous iteration are negligible, we may consider that the algorithm converges practically. In such cases, the iterative updating procedure is stopped, and the parameter estimates obtained are outputted as the final results.

Fig. 1

Diagram of the SAGE algorithm

In the case considered here, we choose the subset of the parameters to be updated in each iteration of the SAGE algorithm to be θ _ℓ =[α _ℓ ,ν _ℓ ,τ _ℓ ], i.e., the parameters of individual paths. The admissible hidden data x _ℓ (f) for estimating θ _ℓ is naturally defined as the contribution of the ℓth propagation path and the noise components, i.e.,

$$ x_{\ell}\left(\,f\right)=\alpha_{\ell} p\left(\,f;\tau_{\ell}, \nu_{\ell}\right)+n'\left(\,f\right). $$

((9))

A SAGE iteration for updating $\hat {\boldsymbol {\theta }}_{\ell }$ consists of two steps, i.e., the so-called Expectation (E-)step and Maximization (M-) step. In the E-step of the ith iteration, the expectation of the loglikelihood of θ _ℓ given y and the estimates of Θ obtained from the ith iteration, denoted with $\hat {\boldsymbol {\Theta }}^{[i]}$, can be calculated as

$${} {\fontsize{9.2pt}{9.6pt}\selectfont{\begin{aligned} \mathrm{E}\left[\Lambda(\boldsymbol{\theta}_{\ell})|\boldsymbol{y}, \hat{\boldsymbol{\Theta}}^{[i]}\right]=&\mathrm{E}\left[-\frac{N}{2}\log 2\pi-\sum\limits_{f=f_{1}}^{f_{N}}\log \mathrm{E}\left[|n'\left(\,f\right)|\right] - \right.\\ &\left.\sum\limits_{f=f_{1}}^{f_{N}} \frac{\left(x_{\ell}\left(\,f\right)-\alpha_{\ell} p\left(\,f;\tau_{\ell}, \nu_{\ell}\right)\right)^{2}}{\mathrm{E}[|n'\left(\,f\right)|^{2}]}\left|\boldsymbol{y}, \hat{\boldsymbol{\Theta}}^{[i]}\right.\right]. \end{aligned}}} $$

((10))

By dropping the constant terms in the right-hand side of (10), it can be shown that

$${} \mathrm{E}\left[\Lambda\left(\boldsymbol{\theta}_{\ell}\right)|\boldsymbol{y}, \hat{\boldsymbol{\Theta}}^{[i]}\right]\propto -\sum\limits_{f=f_{1}}^{f_{N}}\frac{\left(\hat{x}^{[i]}_{\ell}\left(\,f\right)-\alpha_{\ell} p\left(\,f;\tau_{\ell}, \nu_{\ell}\right)\right)^{2}}{\mathrm{E}[|n'\left(\,f\right)|^{2}]}, $$

((11))

where $\hat {x}_{\ell }^{[i]}\left (\,f\right)= \mathrm {E}\left [x_{\ell }\left (\,f\right)|\boldsymbol {y}, \hat {\boldsymbol \Theta }^{[i]}\right ]$ can be calculated as

$$ \begin{aligned} \hat{x}_{\ell}^{[i]}\left(\,f\right) &= y\left(\,f\right) - \sum\limits_{\tiny \begin{matrix}\ell'\neq\ell\\ \!\!\!\ell'\!\!=1\end{matrix}}^{M} \mathrm{E}\left[\alpha_{\ell'} p\left(f;\tau_{\ell'}, \nu_{\ell'}\right)|\hat{\boldsymbol\Theta}^{[i]}\right]\\ &=y\left(\,f\right) - \sum\limits_{\tiny \begin{matrix}\ell'\neq\ell\\ \!\!\!\ell'\!\!=1\end{matrix}}^{M} \hat{\alpha}_{\ell'}^{[i]} p\left(\,f;\hat{\tau}_{\ell'}^{[i]}, \hat{\nu}_{\ell'}^{[i]}\right) \end{aligned} $$

((12))

with $\hat {\alpha }_{\ell '}^{[i]}=\int \alpha _{\ell '} \delta (\alpha _{\ell '} - \hat {\alpha }_{\ell '}^{[i]})\mathrm {d}\alpha _{\ell '}$ and $\hat {\tau }_{\ell '}^{[i]}$, $\hat {\nu }_{\ell '}^{[i]}$ obtained similarly.

For notational convenience, $\boldsymbol {\hat {x}}^{[i]}$ is used in the sequel to represent $\boldsymbol {\hat {x}}^{[i]}=\left [\hat {x}^{[i]}\left (\,f\right); f\in \left [f_{1},\dots,f_{N}\right ]\right ]$.

In the maximization (M-) step of the ith iteration, the estimates $\hat {\nu }_{\ell }^{[i+1]}$, $\hat {\tau }_{\ell }^{[i+1]}$, and $\hat {\alpha }_{\ell }^{[i+1]}$ can be calculated by maximizing the expectation of loglikelihood function obtained in the E-step

$$ \left(\hat{\nu}_{\ell}^{[i+1]}, \hat{\tau}_{\ell}^{[i+1]}, \hat{\alpha}_{\ell}^{[i+1]}\right) = \arg\max\limits_{\nu_{\ell},\tau_{\ell},\alpha_{\ell}}\mathrm{E}\left[\Lambda(\boldsymbol{\theta}_{\ell})|\boldsymbol{y}, \hat{\boldsymbol{\Theta}}^{[i]}\right]. $$

((13))

As shown in Appendix 1, $\hat {\alpha }_{\ell }^{[i+1]}$ can be expressed as a linear function of $\hat {\nu }_{\ell }^{[i+1]}$, $\hat {\tau }_{\ell }^{[i+1]}$ as

$$\begin{array}{@{}rcl@{}} \hat{\alpha}_{\ell}^{[i+1]}=\frac{\boldsymbol{p}\left(\hat{\nu}_{\ell}^{[i+1]}, \hat{\tau}_{\ell}^{[i+1]}\right)^{\mathrm{\!H}}\boldsymbol{W}^{-1}\hat{\boldsymbol{x}}^{[i]}}{\boldsymbol{p}\left(\hat{\nu}_{\ell}^{[i+1]}, \hat{\tau}_{\ell}^{[i+1]}\right)^{\mathrm{\!H}}\boldsymbol{W}^{-1}\boldsymbol{p}\left(\hat{\nu}_{\ell}^{[i+1]}, \hat{\tau}_{\ell}^{[i+1]}\right)} \end{array} $$

((14))

with p(ν,τ)=[p(f;τ,ν);f=f ₁,…,f _N ] being a column vector and W being a diagonal matrix with its diagonal elements equal to E[|n ^′(f)|²],f=f ₁,…,f _N . Inserting (14) to the right-hand side of (13), $\hat {\nu }_{\ell }^{[i+1]}$ and $\hat {\tau }_{\ell }^{[i+1]}$ can be obtained by solving the following maximization problem

$$ \left(\hat{\nu}_{\ell}^{[i+1]}, \hat{\tau}_{\ell}^{[i+1]}\right) = \arg\max\limits_{\nu_{\ell},\tau_{\ell}}\eta\left(\nu_{\ell},\tau_{\ell}\right), $$

((15))

where the objective function η(ν,τ) is shown in Appendix 2 as

$$ \eta(\nu,\tau)=\frac{|\boldsymbol{p}(\nu,\tau)^{\mathrm{H}}\boldsymbol{W}^{-1}\hat{\boldsymbol{x}}^{[i]}|^{2}}{\boldsymbol{p}(\nu,\tau)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\nu,\tau)}. $$

((16))

The amplitude estimate $\hat {\alpha }_{\ell }^{[i+1]}$ can be calculated by using (14).

When the convergence is achieved, e.g., the parameter estimates do not change as the iteration continues, the estimation procedure is suspended, and the parameter estimates obtained in the current iteration are considered to be the final result.

3.1 Discussion of the complexity of algorithm implementation and its influence on the SC device

The complexity of the proposed SAGE algorithm increases along with the number of paths to estimate, the total number of iterations, and the values of B _n , ϱ which determine the number of data samples to be considered when calculating the objective function (16). Reducing the algorithm complexity can be performed in different ways. For example, instead of estimating a large number of multipath components, we may determine an appropriate model order in advance by applying the Akaike Information Criterion [19], and furthermore, the complexity involved in the parameter estimate updating procedure can be reduced by using advanced searching methods [20, 21].

When being implemented in reality, the proposed estimation method requires the SC’s outputs that are the results of filtering the original received sequence by using a LPF of bandwidth B _n , n>1. It can be shown that the relationship B _n ≤f c′ is maintained if n≤γ−1 is selected. Therefore, enlarging the bandwidth of the SC’s output up to B _n with n≤γ−1 would not introduce any additional requirement on increasing the sampling rate at the input of the SC. In such a case, the original complexity in the SC device is not influenced by B _n , 1≤n≤γ−1. The studies in the sequel are conducted under the constraint 1≤n≤γ−1.

4 Simulation study

Simulation studies are conducted to evaluate the performance of the proposed algorithm under the influence of different B _n settings and the fraction ϱ of SC’s output being considered. These two parameters determine how the SC’s output is selected and applied for channel estimation. With larger B _n , more high-frequency components can be involved in parameter estimation, which may improve the estimation accuracy. In the conventional SC-based channel estimation, ϱ=1 is usually adopted. In the proposed algorithm, ϱ<1 can be selected, which allows reducing the observation time required in each snapshot, and consequently, measurements of a time-variant channel can be performed within channel coherence time. The impact of B _n and ϱ on the performance of the proposed algorithm is of importance for understanding the effectiveness of the algorithm. Therefore, we select B _n and ϱ as parameters in the simulation studies.

It is worth mentioning that the conventional channel parameter estimation based on the SC’s output was performed with LPF’s bandwidth set to B ₁=[−f _c /γ,f _c /γ] [7–9]. Therefore, the simulated algorithm’s performance under B=B ₁ can be viewed as the performance achievable when the conventional SC-based estimation is applied. Furthermore, the method proposed here allows estimating the Doppler frequencies of multipath components with SC’s outputs collected in a time period less than that required for calculating a complete estimate of CIR. Since this function is not supported by the conventional SC-based approaches, no references can be provided for comparing with when Doppler frequency estimation results are demonstrated. Without being specifically mentioned, the parameter settings reported in Table 2 are adopted in the simulations. It is worth mentioning that γ=93 is specifically selected here in order to maintain a tractable computational complexity for the simulations. However, γ with larger values may be selected in empirical SC’s applications.

Table 2 Parameter setting of the simulations

The performance of the proposed SAGE algorithm can be investigated from two perspectives, i.e., the RMSEEs in the case where paths are well-separated, and the resolution capability in separating multipath in the case where paths are closely spaced. Sections 4.1, 4.2, and 4.3 are dedicated to the investigation of the RMSEE behavior of the SAGE algorithm in the case with well-separated paths, and Section 4.4 to the resolution capability of the algorithm. When a channel consists of well-separated paths, the received signals of multipath components are mutually orthogonal, and the behavior of the SAGE algorithm can be represented by that of maximum-likelihood estimation (MLE) method in the single-path scenarios [16]. Therefore, the performance of MLE in single-path scenarios are studied in Sections 4.1, 4.2, and 4.3.

4.1 SC’s output with B _n and ϱ as parameters

Figure 2 depicts the magnitude frequency spectrum and the channel power delay profile (PDP) calculated by using of the SC’s output with B _n , n=1 and 7, respectively. The synthetic channel consists of one path with delay of 53 ns, 0 Hz Doppler frequency, and complex amplitude equal to 1·10⁻⁶. The SNR is set to 30 dB. It can be observed from Fig. 2b that for B=B ₁, the channel PDP exhibits a dominant peak located at abscissa close to the true delay. For B=B ₇, i.e., high-frequency components as shown in Fig. 2c are considered, the PDP exhibits a single peak with severe fluctuations on the top. From these results, we can see that the conventional method of delay estimation, which relies on finding the maximum of channel PDP, owns larger estimation errors when the LPF bandwidth increases. In addition, it can be observed from Fig. 2 that in the cases where only a part of the SC’s output is considered, the time-dilated IR obtained with B ₁ can be used for estimating the parameters of a path provided ϱ T>τ ^′, i.e., the acquired part of the SC’s output contains the impulse response of the path. In the case of large B _n , since higher frequency components are considered in the estimation, the part of the IR with ϱ T<τ ^′ may still contain the necessary information for estimating path parameters. It will be shown later than along with the increase of B _n , less ϱ can be considered for estimating path parameters.

Fig. 2

a–d Magnitude frequency spectrum, power delay profile of the CIR calculated by using the low-pass-filtered signals with the bandwidth of B _n , n=1 and 7, respectively. The true path parameters are τ=53 ns, ν=0 Hz, and α=1×10⁻⁶

4.2 Objective functions of delay and Doppler frequency with B _n and ϱ as parameters

4.2.1 η(τ;ν=ν ^′) versus B _n

Figure 3 depicts the objective functions of delay obtained with the LPF bandwidth B _n , n=1,3,5, and 7 considered as a parameter in a single-path scenario. The synthetic channel is generated with the path’s delay equal to 25 ns, Doppler frequency equal to 0 Hz, complex amplitude equal to 1·10⁻⁶. The SNR equals 30 dB. The path’s Doppler frequency ν ^′ is assumed to be known. It can be observed from Fig. 3 that the main-lobes of objective functions exhibit the same zero-to-zero distance in the abscissa regardless of B _n . This indicates that the resolution is maintained the same and irrelevant with B _n . However, the objective function calculated at the correct delay increases along with n. Since the objective function is proportional to the loglikelihood of parameters as shown in Appendix 2, the observation that the maximum of objective function increases along with n implies that enlarging the LPF bandwidth can enhance the likelihood of the delay estimate.

Fig. 3

Objective functions of delay calculated with different B _n in a single-path case with path’s delay equal to 25 ns, Doppler frequency 0 Hz, complex amplitude 1·10⁻⁶, and SNR set to 30 dB. The Doppler frequency is assumed to be known in advance

4.2.2 η(ν;τ=τ ^′) versus B _n

Figure 4 depicts the objective function of Doppler frequency in single-path scenarios. The synthetic channel is generated with path’s complex amplitude equal to 1·10⁻⁶, delay 35 ns, and Doppler frequency 5·10⁴ Hz. The SNR is set to 30 dB. The delay τ ^′ of the path is assumed to be known in advance. Similar observation with that shown in Fig. 3 is obtained, i.e., the objective functions exhibit more peaky main-lobes for larger B _n . This indicates that by enlarging the LPF bandwidth, the likelihood of Doppler frequency can be improved. This is reasonable as more coherent observations are included into the estimation when B _n increases.

Fig. 4

Objective functions of Doppler frequency calculated with different B _n in a single-path scenario where the path’s complex amplitude equals 1·10⁻⁶, delay 35 ns, and Doppler frequency 5·10⁴ Hz. The SNR is set to 30 dB. The path delay is assumed to be known in advance

4.2.3 η(τ;ν=ν ^′) versus ϱ

Figure 5 depicts the objective function η(τ;ν=ν ^′) in the delay domain calculated in a single-path scenario when the length of the SC’s output considered is set with ϱ=1, 1/2, 1/3, and 1/4. The settings for the path and SNR are the same as those adopted in the simulations described in Section 4.2.2. It can be observed from Fig. 5 that the longer the length of SC’s output considered, the higher is the main peak of the objective function. This indicates that the likelihood increases by considering more data at the output of SC. Furthermore, it is observed from Fig. 5 that the zero-to-zero distances of the main peaks of all objective functions are found to be identical, i.e., with the duration of 2 ns, corresponding to the limit of the intrinsic delay resolution set by 1/B with B=500 MHz being the bandwidth of the transmitted signal.

Fig. 5

Objective function of delay calculated with the LPF bandwidth set to B ₆, and ϱ considered as a parameter in a single-path scenario where the path has complex amplitude equal to 1·10⁻⁶, delay 35 ns, and Doppler frequency 5·10⁴ Hz. The SNR is set to 30 dB. The Doppler frequency of the path is assumed to be known in advance in the simulation

4.2.4 η(ν;τ=τ ^′) versus ϱ

Figure 6 depicts the objective function η(ν;τ=τ ^′) of Doppler frequency calculated from the synthetic data in a single-path scenario with SNR of 30 dB. The true delay and Doppler frequency equal 35 ns and 5·10⁴ Hz, respectively. The delay is assumed to be known in the simulations. It can be observed from Fig. 6 that the longer the length of SC’s output considered, the higher is the maximum of the objective function. Meanwhile, the zero-to-zero distance of the main peak of the objective functions is observed to increase along with ϱ. This is reasonable as the total observation span increases when ϱ takes a larger value, resulting in higher resolution in estimating the Doppler frequency. It can be seen that even with a lower resolution, detecting the Doppler frequency is still possible for ϱ<1.

Fig. 6

Objective function of Doppler frequency calculated with the LPF bandwidth set to B ₆, and ϱ considered as a parameter in a single-path scenario where the path has complex amplitude equal to 1·10⁻⁶, delay and Doppler frequency equal to 35 ns and 5·10⁴ Hz, respectively. The SNR is set to 30 dB. The delay is assumed to be known in the simulations

4.2.5 η(τ,ν) versus B _n

Figure 7 illustrates the objective functions η(ν,τ) calculated with the SC’s output in a multi-path scenario with ϱ=1, B _n ,n=1, 7. The synthetic channel consists of a certain number of randomly generated paths with delays chosen from [0,62] ns and Doppler frequencies [−400,400] Hz. The SNR is set to 30 dB with respect to the maximum path power. It can be observed from Fig. 7 that when the B _n increases, the objective function becomes more peaky, and thus, the multipath components can be separated more readily.

Fig. 7

Objective functions η(ν,τ) calculated with a B ₁ and b B ₇ in a multi-path scenario where a certain number of randomly generated paths exist with their delays and Doppler frequencies randomly chosen from [0,62] ns and [−400,400] Hz, respectively. The SNR is set to 30 dB with respect to the maximum path power

Notice that by using the SC’s output that allows generating one CIR, the Doppler frequency estimation resolution is very low due to the short observation duration T. By using the simulation settings in Table 2, the total observation span is calculated to be T=L γ/f _c =31·93·2·10⁻⁹=5.77 μs. Thus, the Doppler frequency estimation resolution is 1/(2T)=43 KHz. Since the differences of Doppler frequencies of paths are empirically much less than 43 KHz, it is important to jointly estimate the delay and Doppler frequency in order to resolve the paths in the delay domain. Furthermore, due to the low Doppler frequency estimation resolution, observations with high SNRs are always preferable in order to obtain less estimation errors. Our simulation results here show that the SNR should be kept beyond 10 and 40 dB in order to obtain RMSEE (ν) less than 10 Hz when the LPF bandwidth B _n is set with n≥5 and n≤3, respectively.

4.3 RMSEEs of delay and Doppler frequency

4.3.1 RMSEE (τ) versus B _n

The benefits by enlarging the LPF bandwidth can also be evaluated by examining the RMSEE of delay and of Doppler frequency. Figure 8 depicts the RMSEE of delay versus the SNR with B _n as a parameter. Single-path scenarios are considered where the Doppler frequencies of the synthetic paths are uniformly distributed within [−1×10⁵,1×10⁵] and the delays of the paths are uniformly selected from 20 to 70 ns. It is worth mentioning that the maximum Doppler frequency for the synthetic paths was set below the maximum frequency that can be measured unambiguously by the system. Totally, 400 Monte-Carlo simulations were conducted. The result obtained by using the conventional estimation method is also illustrated, which searches the estimates of path delays by finding the maxima of channel PDP. This is only feasible when the LPF bandwidth is set to B ₁. Considering that such maxima-searching spectral-based methods are widely used in parameter estimation when the conventional SC is implemented, its performance is taken as a reference of comparison for the parametric estimation algorithm proposed here. It can be observed from Fig. 8 that for a fixed SNR, the RMSEE (τ) when the proposed algorithm is used is at least one order of magnitude lower than those obtained with the conventional method. The worse performance of the conventional method is due to the fact that the time-dilated channel impulse response estimated by using the conventional method consists of 93 samples in the delay domain. In the case where the true path delay is different from integer times of delay samples, by using the maxima-searching method, large estimation errors are resulted. However, in the case where a parametric model-based estimation is performed, searching the maximum of objective function can be performed with refined steps in such a way that more accurate estimates are obtained. In addition, it is observed that RMSEE (τ) decreases when B _n increases, indicating that the estimation accuracy can be improved by taking into account the high-frequency components when the proposed parameter estimation method is applied.

Fig. 8

RMSEE (τ) versus SNR with B _n as a parameter

4.3.2 RMSEE (τ) versus ϱ

As shown in Fig. 5, the resolution capability of the estimator in the delay domain would not be changed when different values of ϱ are considered. This is confirmed by simulations of RMSEE (τ) for two paths with the separation of the two paths in delay taken as a parameter. Figure 9 depicts the results for SNR set to 30 dB and with B ₆. In the simulation, both paths’ delays are randomly selected with specified separation. The fraction ϱ of the IR ranges from $\left [\frac {1}{4}, \frac {1}{3}, \frac {1}{2}, 1\right ]$. Totally 300 snapshots are run for collecting the random samples for computing RMSEEs. It can be observed from Fig. 9 that when the separation Δ τ is larger than 2 ns, the RMSEE (τ) reduces to a stable level, which does not change significantly when Δ τ keeps increasing. This is consistent with the fact that the signal bandwidth of 500 MHz provides the resolution of 2 ns. These results indicate that using the parts of the sliding correlator’s output for parameter estimation would not lead to the reduced resolution. Furthermore, it can also be observed from Fig. 9 that when ϱ increases, the RMSEE decreases. This is due to the reason that the SNR for parameter estimation increases when more observations are included during the calculation of the objective functions.

Fig. 9

RMSEE (τ) of a path 1 and b path 2 versus the separation of two paths in delay with ϱ being a parameter

4.3.3 RMSEE (ν) versus B _n

Figure 10 illustrates the RMSEE of Doppler frequency versus SNR calculated from 250 snapshots in single-path scenarios. The LPF bandwidth B _n changes with n=1,3,5, and 7. The value of ϱ is fixed to 1. It can be observed from Fig. 10 that for fixed SNRs, the RMSEE (ν) decreases as B _n increases. The upper floor of RMSEE (ν) observed for low SNRs in Fig. 10 is due to the Doppler frequency estimation range of [−1,1] KHz specified during the simulations. Furthermore, we observed that the improvements in RMSEE (ν) obtained when the LPF bandwidth changes from B ₁ to B ₃ and from B ₅ to B ₇ are insignificant compared with that resulted when the bandwidth increases from B ₃ to B ₅. This is because the components in the frequency domain do not have the same spectral heights. Although the RMSEE monotonically decreases with respect to the increasing B _n , the decreasing rate is not constant and is actually dependent on the exact range of the abscissa considered.

Fig. 10

RMSEE (ν) versus SNR with B _n as a parameter

4.3.4 RMSEE (ν) versus ϱ

As shown in Fig. 6, the width of main peak of the objective function is enlarged when less fractions of the SC’s output results are taken into account in the estimation. Simulation studies are conducted to verify the resolution ability of the estimator in two-path scenarios when ϱ changes. In the simulations, the bandwidth of the SC output is set to B ₆, and ϱ is set to [1,1/2,1/3,1/4]. The Doppler frequency resolutions κ _ν can be calculated based on κ _ν =1/(2T) to be 7×10⁵, 3.4×10⁵, 2.4×10⁵, and 1.7×10⁵ for ϱ=1/4,1/2,3/4 and 1, respectively. The SNR is set to −10 dB. Figure 11a, b depicts respectively the RMSEE (ν ₁) and RMSEE (ν ₂) for two paths versus their separation in the Doppler frequency domain. The same delay is set for both paths in the simulations and is assumed to be known in advance in parameter estimation. It can be observed from Fig. 11a, b that both RMSEE (ν ₁) and RMSEE (ν ₂) decrease when Δ ν increases. Practically, we can define the empirical resolution as the separation of two paths beyond which the resultant RMSEE (ν) for both paths becomes stabilized. It can be observed from Fig. 11 that the empirical resolutions are consistent with the theoretical intrinsic resolutions which are defined to be the inverse of the total observation span, and moreover, the empirical resolutions are found to increase when ϱ decreases. These results show that it is possible to estimate Doppler frequency using parts of the SC’s output, and that the length of the SC’s output determines the resolution of separating paths in the Doppler frequency domain. It is worth mentioning that significant fluctuations can be observed for RMSEE (ν ₁) and RMSEE (ν ₂) graphs when the Doppler frequency separation is less than the intrinsic resolution of the estimator. This is due to the biases in the Doppler frequency estimates. When Δ ν between two synthetic paths is less than the resolution, the maximum of the objective function calculated for estimating the Doppler frequency of the first path in the initialization step is usually found between the true Doppler frequencies. Although the SAGE algorithm can change the estimates with more iterations, such biases may still exist in the final estimation results. In addition, since the biases have different values depending on Δ ν, the RMSEE graphs obtained with Δ ν less than the intrinsic resolution exhibit significant fluctuations.

Fig. 11

RMSEE (ν) of a path 1 and b path 2 versus the separation of two paths in Doppler frequency with ϱ being a parameter

4.3.5 RMSEE (τ) versus B _n and ϱ

The aforementioned investigations focus on the behavior of the estimator with respect to either B _n or ϱ. We now try to evaluate how ϱ and B _n jointly influence the performance of the proposed estimation algorithm. Figure 12 illustrates RMSEE (τ), i.e., the RMSEE of delay in a single-path scenario with B _n varying within the range n=1,2,3, and the fraction ϱ of the SC’s output taking values among $\left [\frac {1}{4}, \frac {1}{3}, \frac {1}{2}\right ]$. The SNR of 10 dB is considered in the simulations. Totally, 300 snapshots are applied to get the RMSEE (τ) graphs. The true parameters of the path are τ=25 ns, ν=0 Hz, |α|=3. The estimation of delay is perfromed under the assumption that the Doppler frequency is known in advance. It can be observed from Fig. 12 that when we take only a part of the SC’s output for estimating the parameters, only for the bandwidth $B\geq n\frac {f_{c}}{\gamma }$ with n=3, the parameters can be estimated.

Fig. 12

RMSEE (τ) calculated in single-path scenarios with known Doppler frequency

4.3.6 RMSEE (ν) versus B _n and ϱ

Figure 13 depicts RMSEE (ν), i.e., the RMSEE of Doppler frequency with the bandwidth B _n being variable within the range $\left [n\frac {f_{c}}{\gamma }; n= 1,2,3\right ]$, and the fraction ϱ of the SC’s output taking values among $\left [\frac {1}{4}, \frac {1}{3}, \frac {1}{2}\right ]$. Simulation settings for the path parameters as in the simulations for generating Fig. 12 are applied. The estimation of Doppler frequency is conducted under the assumption that the delay is known in advance. It can be observed from Fig. 13 that when we take only a part of the SC’s output for estimating the parameters, for the bandwidth n≥3, the parameter estimation returns RMSEEs that start to be stabilized. This is consistent with the observations from Fig. 12. It can be also observed from Fig. 13 that the RMSEE (ν) graphs fluctuation when n varies. A possible reason for this effect is that a larger n does not necessarily lead to more signal components included into the estimation. From Fig. 2, we observed that when n increases from 1 to 7, the number of the mainlobes of signal components increases from 1 to 5, which implies that the signal contribution to the observations applied in parameter estimation does not increase linearly with respect to n. We postulate that this uneven increments of signal components generates the fluctuations of the RMSEE (ν) shown in Fig. 13.

Fig. 13

RMSEE (ν) calculated in single-path scenarios with known delay and SNR equal to 10 dB

4.4 The SAGE performance in two-path scenarios

4.4.1 RMSEE (τ) and RMSEE (ν) versus B _n

The performance of the derived SAGE algorithm is evaluated in a two-path scenario, where the parameters of the paths are (τ ₁,ν ₁,α ₁)=(22 ns,−40 Hz,3) and (τ ₂,ν ₂,α ₂)=(28 ns,40 Hz,1), respectively. The noise components are added with N ₀=max(|α ₁|²,|α ₂|²)10^−ζ/10 where max(a,b) returns the maximum of the given arguments a and b, ζ denotes the SNR in dB. To limit the simulation times, the SAGE algorithm was executed to estimate the parameters of two paths within maximum 5 iterations. Since the true paths are set with different magnitudes, the path estimated by the SAGE algorithm with larger magnitude is considered to the estimate of the first path, and the other estimated path is the estimate of the second path, i.e., the weaker path.

Figure 14 depicts the RMSEEs of delay and Doppler frequency for the two paths versus SNR obtained from 250 simulation snapshots. It can be observed from Fig. 14 that the RMSEEs of parameters for the second path are always larger than their counterparts for the first path. This is reasonable since the path 2 has lower SNR than the path 1 when synthetic channels were generated. Furthermore, by comparing Fig. 14a, b, we observe that the decrease of the RMSEE of Doppler frequency for the second path when the LPF bandwidth increases from B ₅ to B ₇ is more significant than that obtained for the first path. We postulate that this effect of larger improvement in parameter estimation obtained for weaker paths by increasing the LPF bandwidth is due to two reasons, i.e., a higher objective of the parameters can be resulted when B _n is enlarged, and additionally, the interference cancelation in the E-steps of the SAGE algorithm can be performed more efficiently, especially for paths with lower power.

Fig. 14

a–dRMSEEs of delay and Doppler frequency in two-path scenarios

4.4.2 RMSEE (τ) and RMSEE (ν) versus ϱ

The performance of the derived SAGE parametric estimation algorithm developed is investigated in two-path scenarios. In the simulations, the bandwidth B _n is equal to $6\frac {f_{c}}{\gamma }$, and the fraction ϱ is set to the values within $\left [\frac {1}{4}, \frac {1}{3}, \frac {1}{2}\right ]$. The SNR varies from 0 to 30 dB in step of 10 dB in the simulations. The true path parameters are set to [τ ₁,ν ₁,α ₁]=[22 ns,−40 Hz,2] and [τ ₁,ν ₁,α ₁]=[28 ns,40 Hz,1]. The RMSEE graphs were generated by using 300 snapshots. It can be observed from Fig. 15 that when the fraction ϱ increases, the RMSEEs reduce. Since less ϱ may also influence the resolution in Doppler frequency domain, the improvement in RMSEE(ν) attributed to a larger value of ϱ is more significant than that observed for RMSEE(τ).

Fig. 15

a–d RMSEEs of delay and Doppler frequency for two paths in two-path scenarios with B ₆ and ϱ considered as a parameter

4.4.3 Performance of the SAGE algorithm in a multipath propagation scenario

The performance of the SAGE algorithm is investigated in the case where a channel consists of 10 paths. The average SNR equals 10 dB. The SAGE algorithm was set to estimate 10 paths from the received signals. The results show that running the SAGE algorithm for 10 iterations is sufficient for observing a stabilized likelihood of parameter estimates. Figure 16 depicts a comparison of synthetic and estimated multipath components obtained after 10 iterations. It is obvious by comparing Fig. 16a, b that the parameters of the estimated paths are not exactly the same as those of true paths. This is reasonable since the SAGE algorithm, which approximates the MLE with iterative procedures, has an inherent limitation that the estimates may lead to a local maximum of likelihood that is usually determined by the initialized parameter estimation instead of a global maximum. Furthermore, the erroneous results may also attribute to the limits of the intrinsic resolutions in delay and Doppler frequency domains, since paths closely spaced by distances less than the resolutions cannot be resolved accurately by the SAGE algorithm. Furthermore, we can also observe that paths estimated with significant magnitudes in Fig. 16b appear in the vicinity of their counterparts observed in Fig. 16a. Figure 17 illustrates the normalized delay-Doppler frequency power spectra (PS’s) P _Bartlett(τ,ν) calculated by using the Bartlett beamforming technique [17] of the received signals, $\hat {P}_{\text {Bartlett}}(\tau,\nu)$ of the reconstructed signal calculated based on the SAGE estimation results, and $\tilde {P}_{\text {Bartlett}}(\tau,\nu)$ of the residual signals calculated by subtracting the reconstructed signals from the original received signals. For comparison convenience, these PS’s are all normalized by the maximum of P _Bartlett(τ,ν) and represented in dB in Fig. 17. It can be observed from Fig. 17 that $\hat {P}_{\text {Bartlett}}(\tau,\nu)$ is consistent with P _Bartlett(τ,ν) especially in the portions of larger spectral height, and the maximum of $\tilde {P}_{\text {Bartlett}}(\tau,\nu)$ for residual signals is 34 dB below the maximum of P _Bartlett(τ,ν). These results demonstrate that the proposed SAGE algorithm is capable of extracting dominant components in the multipath channel, although the estimated paths may not have exactly the same parameters as the true paths due to the existence of noises, the limited resolutions caused by finite signal bandwidth and observation spans, as well as the inherent limitation of the SAGE algorithm.

Fig. 16

Comparison of the a synthetic and bestimated multipath components

Fig. 17

Comparison of the normalized delay-Doppler frequency power spectra a P _Bartlett(τ,ν) of the received signals, b $\hat {P}_{\text {Bartlett}}(\tau,\nu)$ of the reconstructed signal calculated based on the SAGE estimation results, and c $\tilde {P}_{\text {Bartlett}}(\tau,\nu)$ of the residual signals calculated by subtracting the reconstructed signals from the original received signals

5 Conclusions

In this contribution, a parametric generic model was proposed to describe the output of the sliding correlator (SC) which is usually utilized to calculate the time-dilated wideband propagation channel impulse responses (CIRs). Based on the model proposed, a Space-Alternating Generalized Expectation-maximization (SAGE) algorithm was derived for extracting the delays, Doppler frequencies, and complex attenuations of multipath from the SC’s output that contains only one observation of time-dilated CIR. Simulation results have shown that the conventional constraint that only the low-frequency component in the SC’s output is applicable for channel estimation is unnecessary when the proposed estimation method is used. Furthermore, more high-frequency components considered, the higher the estimation accuracy can be achieved. Compared with the conventional approach which estimates the channel based on the time-dilated CIR, the proposed method is applicable not only for estimating the multipath’s Doppler frequencies but also returns more accurate estimates than the conventional method, e.g., the delay estimation errors are at least one order of magnitude less than those obtained by using the conventional method. Simulation results also demonstrated that the root mean squared estimation errors (RMSEEs) can be reduced by enlarging the bandwidth of a low-pass-filter (LPF) applied in the SC. When only a part of the SC’s output is available, the parameters can still be estimated provided the bandwidth of the LPF is no less than three times of the transmitted signal bandwidth divided by the sliding factor. In cases where only fractions of SC’s output is considered for estimation, the RMSEEs increase along with the data amount due to the improved output signal to noise ratio and the enhanced resolution capability particularly in the Doppler frequency domain. These results revealed the potential of applying the proposed high-resolution method in the SC-based parameter estimation for mm-wave wideband channel characterization.

6 Appendix 1: Derivation of (14)

As shown in the right-hand side of (11), the original loglikelihood function by dropping constant terms can be rewritten as

$$\begin{array}{@{}rcl@{}} L(\boldsymbol{\theta}_{\ell}) &=& -\sum\limits_{f=\,f_{1}}^{f_{N}}\frac{\left(\hat{x}^{[i]}_{\ell}\left(\,f\right)-\alpha_{\ell} p\left(\,f;\tau_{\ell}, \nu_{\ell}\right)\right)^{2}}{\mathrm{E}\left[|n'\left(\,f\right)|^{2}\right]} \end{array} $$

((17))

$$\begin{array}{@{}rcl@{}} &=& -\left(\hat{\boldsymbol{x}}_{\ell}^{[i]} - \alpha_{\ell} \boldsymbol{p}(\tau_{\ell},\nu_{\ell})\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\left(\hat{\boldsymbol{x}}_{\ell}^{[i]} - \alpha_{\ell} \boldsymbol{p}(\tau_{\ell},\nu_{\ell})\right),\\ \end{array} $$

((18))

where W is a N×N diagonal matrix calculated as

$${} {\fontsize{7.9pt}{9.6pt}\selectfont{\begin{aligned} \boldsymbol{W} \,=\, \left[\!\begin{array}{ccccc} \mathrm{E}\left[|n'\left(\,f_{1}\right)|^{2}\right] & 0 & 0 & \dots & 0 \\ 0 & \mathrm{E}\left[|n'\left(\,f_{2}\right)|^{2}\right] & 0 & \dots & 0\\ 0 & 0 & \mathrm{E}\left[|n'\left(\,f_{3}\right)|^{2}\right] & \dots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \dots & \mathrm{E}\left[|n'\left(\,f_{N}\right)|^{2}\right]\end{array}\!\right] \ . \end{aligned}}} $$

((19))

When the parameters τ _ℓ , ν _ℓ are given, α _ℓ can be determined by solving the following equation

$$ \frac{\partial L(\boldsymbol{\theta}_{\ell})}{\partial \alpha_{\ell}}=0. $$

((20))

From (18) it is easy to show that $\frac {\partial L(\boldsymbol {\theta }_{\ell })}{\partial \alpha _{\ell }}$ can be calculated as

$${} {\fontsize{8.2pt}{9.6pt}\selectfont{\begin{aligned} \frac{\partial L(\boldsymbol{\theta}_{\ell})}{\partial \alpha_{\ell}} &= \frac{\partial\left(\hat{\boldsymbol{x}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1} \hat{\boldsymbol{x}}_{\ell}^{[i]}}{\partial \alpha_{\ell}} - \frac{\partial\alpha_{\ell}\left(\hat{\boldsymbol{x}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) }{\partial \alpha_{\ell}} \\ &- \frac{\partial\alpha_{\ell}^{*}(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\hat{\boldsymbol{x}}_{\ell}^{[i]}}{\partial \alpha_{\ell}} + \frac{\partial\alpha_{\ell}^{*}(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})\alpha_{\ell}}{\partial \alpha_{\ell}}\\ &=-\left(\hat{\boldsymbol{x}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) + \hat{\alpha}_{\ell}^{*}(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \end{aligned}}} $$

((21))

Applying (21) in (20) yields

$${} \left(\hat{\boldsymbol{x}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) - \hat{\alpha}_{\ell}^{*}(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})=0, $$

((22))

which further leads to the expression of $\hat {\alpha }_{\ell }^{*}$ as a function of τ _ℓ , ν _ℓ , and $\hat {\boldsymbol {x}}_{\ell }^{[i]}$:

$$ \hat{\alpha}_{\ell}^{*} = \frac{\left(\hat{\boldsymbol{x}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}{ (\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}\left(\tau_{\ell},\nu_{\ell}\right)}. $$

((23))

By taking the complex conjugates of both sides in (23), (14) is finally obtained.

7 Appendix 2: Derivation of (16)

Substituting α _ℓ in (21) by (14) yields for L(θ _ℓ )

$$\begin{array}{@{}rcl@{}} L(\boldsymbol{\theta}_{\ell}) &=&\! \left(\!{\boldsymbol{\hat{x}}}_{\ell}^{[i]}\!\right)^{\mathrm{H}}\!\boldsymbol{W}^{-1} {\boldsymbol{\hat{x}}}_{\ell}^{[i]} - \frac{\left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}{\boldsymbol{\hat{x}}}_{\ell}^{[i]}({\boldsymbol{\hat{x}}}_{\ell}^{[i]})^{\mathrm{H}} \boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) }{ \left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}\\ &-&\!\frac{\left({\boldsymbol{\hat{x}}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}{\boldsymbol{\hat{x}}}_{\ell}^{[i]}} {\left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell})\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})} \\ &&\quad+\frac{\left({\boldsymbol{\hat{x}}}_{\ell}^{[i]}\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}{ \left(\boldsymbol{p}\left(\tau_{\ell},\nu_{\ell}\right) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})} \left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})\\ &&\qquad\quad\frac{(\boldsymbol{p}\left(\tau_{\ell},\nu_{\ell}\right))^{\mathrm{H}}\boldsymbol{W}^{-1}{\boldsymbol{\hat{x}}}_{\ell}^{[i]}} {(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}. \end{array} $$

((24))

It is easy to show that the last two terms in (24) are identical to each other with opposite signs. Thus, (24) can be rewritten as

$${} {\fontsize{8.4pt}{9.6pt}\selectfont{\begin{aligned} L(\boldsymbol{\theta}_{\ell}) \,=\, \left(\!{\boldsymbol{\hat{x}}}_{\ell}^{[i]}\!\right)^{\mathrm{H}}\!\boldsymbol{W}^{-1} {\boldsymbol{\hat{x}}}_{\ell}^{[i]} - \frac{\left(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \right)^{\mathrm{H}}\boldsymbol{W}^{-1}{\boldsymbol{\hat{x}}}_{\ell}^{[i]}\left(\!{\boldsymbol{\hat{x}}}_{\ell}^{[i]}\!\right)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) }{ \bigl(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}) \bigr)^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}. \end{aligned}}} $$

((25))

The term $\left ({\boldsymbol {\hat {x}}}_{\ell }^{[i]}\right)^{\mathrm {H}}\boldsymbol {W}^{-1} {\boldsymbol {\hat {x}}}_{\ell }^{[i]}$ in the right-hand side of (25) is constant with respect to θ _ℓ . By dropping this constant term, we obtain

$$ L(\boldsymbol{\theta}_{\ell}) \propto - \frac{\bigl|(\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}{\boldsymbol{\hat{x}}}_{\ell}^{[i]}\bigr|^{2}}{ (\boldsymbol{p}(\tau_{\ell},\nu_{\ell}))^{\mathrm{H}}\boldsymbol{W}^{-1}\boldsymbol{p}(\tau_{\ell},\nu_{\ell})}. $$

((26))

From (26) it is obvious that minimization of L(θ _ℓ ) with respect to θ _ℓ is equivalent with maximization of an objective function η(τ _ℓ ,ν _ℓ ) defined as shown in (16).