Methods for mutual time delay estimation of wideband signals based on nonlinear digital filtering

. Methods for mutual time delay estimation of wideband signals propagating in satellite communication systems are proposed. The signals are propagated in different channels and received with low signal-to-noise ratio. A characteristic feature of satellite channel is the presence of the Doppler Effect, which leads to a shift and scaling the signal spectrums. The proposed approaches are based on the separation of narrow-band channels from the studied signals, using non-linear digital filtering algorithms in each channel, and subsequent optimal (correlation) processing. The accuracy of the proposed methods and the reliability of the determination of time delay are investigated.


Introduction
One of the most important problems in modern satellite communication systems is the problem of determining position of an emitting object in real time by passive direction finding techniques. For example, the task is relevant for switching and establishing a communication session between spacecraft and receiver located near the Earth ground. To establish a communication session it is necessary to form a radiation pattern using a phased antenna array located on spacecraft. To set direction of beam it is necessary to know location of emitting object on the Earth and determine it in real time. Another example of passive direction finding methods application is the problem of designing satellite search and rescue systems.
One of the most common methods for determining location of radiation source is the range-difference method. This method requires an estimate of mutual time delays of distorted copies of emitter signal [1]. Distorted copies of signal are result from the propagation in different channels with a low signal-to-noise ratio (SNR).
Wideband signals are widely used in modern satellite communication systems due to high reliability of the information transmission and robustness to noise. In particular signals with pseudorandom frequency hopping (Frequency hopping spread spectrum technology, FHSS) and signals with orthogonal frequency-division multiplexing (OFDM) are used.
Currently a sufficient amount of algorithms for estimating the mutual time delay has been developed. The basis of these algorithms is the maximum likelihood method, which can give robust estimate. However, direct use of this approach is ineffective when working with wideband signals in presence of the Doppler Effect. By effectiveness of the algorithm we mean the ability of the method to give consistent real time estimates in low SNR conditions and presence of the Doppler Effect with wide range of a priory uncertainty of mutual time delays. The effectiveness of most of the algorithms developed to date is significantly reduced when there is spectrum scaling which is a consequence of the Doppler Effect influence on wideband signal.
Therefore, the problem of developing methods to obtain a reliable estimate of the mutual time delay of wideband signals in real time is actual. In this paper methods for determining the mutual time delay based on nonlinear filtering are developed and investigated.
2 Methods for estimating the mutual time delay of wideband signals

The problem of time delay estimation
Let us assume that two signals s(t) and () st propagate in different channels. The propagation model can be written as follows: The signal x(t) is generally unknown. Signal x(t) is received in one of the channels. Signal () st is delayed by time τ and represents a distorted copy of the signal s(t). Signals ξ(t) and η(t) are noises, which are uncorrelated with the signals. Signal s(t) is taken as a reference, since its signal-to-noise ratio is sufficiently large (about +10 dB). Typically, the reference signal is the signal received from a geostationary satellite. Another signal () st is taken as observed signal, and its signal-to-noise ratio could be less than zero. The task is to estimate the mutual time delay τ.
If the propagating signals are the realizations of stationary Gaussian processes with zero mean, the asymptotically optimal method in the sense of maximum likelihood is calculation of the cross-correlation between pre-filtered signals 1 () st and 2 () st ("generalized crosscorrelator") [2]: where T is the duration of the reference signal s 1 (t). Pre-filtering is carried out on the basis of a priori information about noise spectrum (if it is available) in order to increase the SNR. The mutual time delay is determined by the position the main peak of the cross-correlation: However, in the communication systems with mobile objects, particularly used in space communication, the correlation methods require a compensation for the shift of the signal spectrums caused by the influence of the Doppler Effect. The classical method of time delay estimation that does not require any additional knowledge about the signal structure is based on the calculation of the cross-ambiguity function [3]: For narrowband signals, the time delay and Doppler shift can be estimated from the position of the main peak of the absolute value of the cross-ambiguity function in timefrequency dimension. However, it is important that the Doppler Effect leads to two kind of distortion of signal spectrums. First, it is the shift of the carrier frequency f 0 by the constant value of Δf ≈ α•f 0 , where α is a factor depending on relative velocity of emitter and receiver. The second effect is scaling the signal spectrum by factor (1+α). Since the characteristic value of α for the satellite system (assume that the satellite moves in a high elliptical orbit) is in order α ~ 10 -5 , for the signals with bandwidth about B ~ 100 MHz the value of spectrum scaling is about 1 kHz. This value is not negligible compared the frequency shift (in order Δf ≈ 10 kHz for f 0 ~ 1 GHz).
For wideband signals, the cross-ambiguity function cannot accurately compensate for the frequency shift due to the strong influence of the effect of scaling the spectrum. The degree of intensity of the main peak of the ambiguity function of wideband signals decreases sharply with the increase in a value of the Doppler shift of the carrier frequency and scale of the spectrum. Consequently, the direct calculation of the ambiguity function of the ultra-wideband signals does not allow for a consistent estimate of the mutual time delay in low SNR conditions [4].

Modified algorithm for estimating the time delay
To determine the mutual time delay of wideband signals we propose a method based on separation of M narrowband channels from received wideband signals using a set of bandpass filters [4,5]. Consequently, digital nonlinear filtering is applied to signals in these channels. The general scheme of the algorithm is presented in Fig. 1. s(t) and () st are input signals between which it is necessary to determine the delay. H i (f) is a gain characteristic of the band-pass filter to separate i-th narrowband channel. G i is nonlinear filter, R i is the cross-correlation function of i-th pair of outputs of the nonlinear filter. R is the resulting time distribution (cross-correlation) obtained as a result of addition of the R i functions over all narrowband channels, i = 1..M.  The bandwidth and central frequency of each channel are set at the initial stage of the algorithm and are the parameters of band-pass filters. Increasing the computational efficiency of the method and saving the required memory can be achieved by choosing the optimal number M of channels.
Since signals in channels are narrowband, the algorithm of optimal processing can be applied. At the same time the spectrum scaling in channels can be ignored and it is necessary to compensate the shift of central frequency of signal. To compensate the shift algorithm for calculating the cross-ambiguity function (4) can be applied. However, this algorithm has high computational complexity due to wide range busting of time and frequency shifts.
The methods based on nonlinear filter avoid compensation for unknown frequency shift of signals in narrowband channels. Nonlinear filtering is a method based on replacing the samples of the original signals with samples of some function depending on the "instantaneous" frequency of the signals followed by optimal correlation processing of the obtained functions.
Addition of the cross-correlation functions over all narrowband channels is applied to increase the intensity of the main peak of the resulting time distribution. The crosscorrelation R i has the main peak with low intensity due to low SNR, short information part in channel (if the signals are i.e. FHSS signals) [6], low information transfer rate in channel and influence of adjacent channels (if the signals are i.e. OFDM signals) [4].

Nonlinear filter based on Pisarenko harmonic decomposition
This method is similar to signal demodulation procedure and reduces to processing the original signals with an adaptive digital filter. This filter replaces signal samples with samples of another function that we call the "current frequency function". To form this function a "sliding" calculation of the autocorrelation function and Pisarenko harmonic decomposition (PHD) are used.
To find "instantaneous" frequencies by PHD it is necessary to compile the Toeplitz correlation matrix R xx and to solve the matrix equation of the following form [7]: where σ n 2 is the noise variance, A is the parameter vector of the autoregressive model which is an eigenvector of the R xx matrix at the same time A=(a 0 , a 1 ,…,a p ), p is an order of the autoregressive model. It is known that the noise variance corresponds to the minimal eigenvalue of the matrix R xx if the signal is the superposition of sinusoids with Gaussian white noise.
The eigenvector A=(a 0 , a 1 ,…,a p ) can be determined from eq. (5) by using singular value decomposition (SVD) method. From the coefficients of this vector the polynomial on z are formed: The complex roots of the polynomial (6) determine the current frequency [7]: where T s is the sampling period.
Since the signal in narrowband channel can be a harmonic modulated signal let us assume that the order p = 3. This is a pretty rough assumption therefore the main peak of the cross-correlation function in one channel has low intensity. However, the proposed algorithm based on obtaining the resulting time distribution by addition of crosscorrelations over all channels allows us to get a well-defined main peak and blurred side maxima.
In case p = 3 the polynomial (6) is the quadratic equation for z in complex numbers: moreover A=(1, a 1 /a 0 , a 2 /a 0 ). The "instantaneous" frequency value of signal is determined according to the following expression: So, the process of forming the "current frequency function" can be divided into several stages: 1. The length of "sliding" window L is selected (the number of signal samples from which the autocorrelation function is calculated). 2. A Toeplitz 3x3 matrix R xx is compiled from the samples of signal autocorrelation function in "sliding" window. 3. The minimum eigenvalue λ min (corresponding to the noise variance) and the corresponding eigenvector A=(1, a 1 ', a 2 ') are determined (i.e. by SVD algorithm). 4. The quadratic equation (8) is made up of coefficients of vector A and the value of frequency is found by (9). 5. The window moves by 1 sample, and go to step 2. The described algorithm calculates the "current frequency function" for each narrowband channel. So we obtain two sets of the "current frequency function" corresponding to signals s(t) and () st . The procedure of nonlinear filtering qualitatively resembles demodulation process and allows us to detect frequency deviations from constant frequency. Moreover, these deviations are associated with phase gaps in signals (information bit change). Therefore, the cross-correlator (2) can be applied for the "current frequency function" and for increasing the intensity of the main peak the resulting time distribution is obtained by summing the cross-correlations over all channels. The time delay between signals s(t) and () st is estimated by the main peak position of the resulting time distribution (3).

Quadratic filter based on Capon approach
A quadratic filter was proposed in paper [8]. This filter based on a generalization of the Capon minimum variance approach. The Capon approach consists in constructing a linear filter minimizing output signal variance when transmission coefficient on given frequency equals one: where the index H denotes Hermitian conjugation, c is the vector of filter coefficients, e(f 0 ) is the vector of complex exponentials whose components are equal e k (f 0 )=exp(j2πf 0 k), k is the index of component, R xx is the autocorrelation matrix of harmonic signal with frequency f 0 . An analytical solution to this problem is obtained by conditionally minimizing the variance of the filter output by the Lagrange multiplier method. The solution can be written as follows [7,8]: where R xx -1 is the inverse matrix to the autocorrelation matrix R xx . It is assumed that the autocorrelation matrix is non-degenerate and the number L of filter coefficients is determined by the rank of this matrix. Increasing of L when calculating the autocorrelation matrix from the processed signal makes the estimation of the Capon filter coefficients unstable. It creates obstacles to the practical implementation of the minimum dispersion method in signal filtering problems. Nevertheless, the idea of the Capon method can be serve as the basis for constructing nonlinear filters devoid of the limitations noted above.
A quadratic filter based on the linear filter (11) can be constructed. The output signal of the filter is determined by the expression [8]: where the exponential factor is responsible for taking into account the initial phase of processed signal. In practice accounting for this factor is quite difficult. Failure to do so leads to a slight (for low SNR) deterioration in the output signal form. The calculation of the inverse matrix for an arbitrary order of the filter should be replaced by the calculation of the pseudoinverse Moore-Penrose matrix. The denominator of expression (12) is a normalization factor and can be omitted for ease of calculation. As a result the filter output signal is obtained by the following expression [8]: where R xx # (f 0 ) is the Moore-Penrose matrix of the autocorrelation matrix of the harmonic signal on frequency f 0 .
The quadratic filter (13) is the narrowband filtering system so it can be the basis for developing algorithms for narrowband signals processing.
Similar to section 2.2.1 the output signal of quadratic filter can be considered as "current variance function". The procedure of calculation this function also qualitatively resembles the demodulation process. Therefore, the cross-correlator (2) can be applied for the "current variance function" and the time delay is estimated by the position of the main peak in resulting time distribution calculating as a sum of cross-correlation functions over all narrowband channels.

Simulation
Computer simulations are carried out a study characteristic of the proposed algorithms for time delay estimation. OFDM-signals were created with following parameters:  bandwidth B=400 MHz;  carrier frequency f 0 =1 GHz;  number of subcarriers N=512. When simulating the Doppler shift of carrier (central) frequency Δf 0 =10..50 kHz was used, and the coefficient of scaling the band α=Δf/f 0 was calculated. Modelling of scaling the spectral band of signal was carried out by performing oversampling of the original signal at (1+α) times based on spline interpolation. Next, the proposed algorithms are applied for created signals. For this, signals were separated into M narrowband channels via linear band-pass filtering. The signal bandwidth in one channel was chosen equal to B k =5 MHz, so scaling the spectrum in channel is negligible compared to shift of the spectrum.
A study was conducted to determine the optimal number of channels. To determine the degree of reliability of the time delay estimate we propose to use a dimensionless and energy-independent criterion: where R(Δt) is the resulting time distribution obtained by the proposed algorithm (Fig. 1), max|R(Δt)| is the maximal absolute value of this distribution, () Rt  is the mean absolute value and var|R(Δt)| is the variance of the distribution.
The criterion also determines intensity of the main peak in resulting time distribution. When the time distribution has well-defined main peak, the value of criterion is large. On the other hand when the level of side maxima is large, the main peak has low intensity and the criterion value is small because of the large value of the variance. Figure 2 shows the relation between the criterion C and the number of narrowband channels at various values of SNR for the algorithm based on the quadratic filter. As could be seen from Fig. 2 the value of criterion increases with increasing the number of narrowband channels M. When M≥15 the function C(M) reaches saturation. Based in this, for each channel bandwidth we can determine the effective value of the number M. So, for the B k =5 MHz let us assume than it is sufficient to separate 15 narrowband channels to obtain a well-defined main peak in the resulting time distribution. For the algorithm based on PHD a similar dependency was obtained ( Figure 3). As could be seen from Fig. 3 the value of M=12 is enough to obtain a well-defined main peak in the resulting time distribution.
For the algorithm based on quadratic filtering it is necessary to determine the optimal value of the filter length, since its increase leads to a decrease in the computational efficiency of the algorithm and a decrease in length leads to worsening of the main peak intensity in the resulting time distribution. The obtained dependency C(L) is presented in Figure 4.  Analyzing the data presented in Fig. 4 it can be seen that the dependence can be divided into two parts -increasing and decreasing. The maximum value of criterion is at a filter length value of L=7. This maximum can be explained by the fact that a certain filter length is necessary for the algorithm to work optimally, because a filter that is too short is built on a small amount of information about the signal, and too long introduces instability into the algorithm work. Thus, for further estimation of the mutual time delay using this algorithm the filter length L=7 is used.
For testing, the stability of the algorithms on low SNR a statistical research was performed. A series of experiments was carried out to evaluate the mutual time delay at various SNR values. The experiment consisted of calculating the time delay and fixing the value of estimate falling within the confidence interval and compliance with the Neumann-Pearson criterion (with the probability of false alarm 0.01). The length of one information symbol in input signal is selected as a confidence interval. An average over 1000 experiments for each value of SNR. The probability P of falling the confidence interval is calculated as the ration of successful experiments to the total number of experiments. Figure 5 shows the obtained dependencies P(SNR) for the algorithm base on quadratic filter, the algorithm based on PHD, and (for comparison) the algorithm based on calculating the cross-ambiguity function for signals in narrowband channels. The algorithms based on calculation the cross-ambiguity function allows us to determine the mutual time delay at SNR up to -16 dB. However, the disadvantage is the low computational efficiency and the requirement of a large a number of resources due to the complexity of calculating the cross-ambiguity function. This implies the inability to use this algorithm in real time. The proposed algorithms also allow us to reliably determine the mutual time delay at low SNR (up to -12 and -13 dB for the algorithm based on quadratic filter and the algorithm based on PHD respectively), but at the same time they have a high speed of operation and do not require a large amount of computing resources. Comparing the dependencies for the algorithm base on quadratic filtering and the algorithm based on PHD we can see that the latter works more stably.

Conclusion
The efficient algorithms for time delay estimation of ultra-wideband signals are proposed. The algorithm are based on splitting the signals into a set of narrowband channels and the subsequent nonlinear filtering. A filter based on Pisarenko harmonic decomposition and quadratic filter based on Capon approach are considered. Summing cross-correlation functions of the nonlinear filter outputs across all channels can significantly increase the probability of reliable time delay determination.
The proposed algorithms are stable at low SNR and they are more computation effective than the algorithm based on calculation the cross-ambiguity function.
The proposed algorithms can be used as a basis of the algorithm for determining the location of the radiation source by the range-difference method.