Investigation of the time delay difference estimator for FMCW signals Mikhail V. Ronkin and Aleksey A. Kalmykov Ural Federal University, Yekaterinburg, Russia MVRonkin@gmail.com Abstract. The paper deals with the problem of development of the effective time delay difference estimator of beat signals, obtained by pro- cessing the frequency modulated signals. Such beat signals can be formed by the heterodyne scheme of receiver of linear frequency modulated sig- nals with internal coherence. The described problem is actual in many radar applications such as radio direction finding. The paper is carried out the analysis of Cramer-Rao lower bound for delay as signal param- eter and its comparison to frequency and initial phase. On the base of the considered approach the estimator of delay difference of beat signals is proposed in the class of smooth parabolic finite difference estimators. Investigation of the proposed method is carried out. Advantages of the developed estimator are shown in comparison with other techniques in such class in the areas of variance and influence of the parasitic signals on the bias. Keywords: signal processing, FMCW signals, chirp signals, beat sig- nals, delay measurements, radars, delay estimators 1 Introduction The problem of time delay difference estimation of received signals is actual in many applications of radar, sonar systems and also in ultrasonic, laser and other measurement tasks [1], [2], [3], [4], [5]. The main problem formulation of time delay difference measurement is estimation of ∆τ from the resulted samples of two received signals [5] s1 (t) = h1 (t) ∗ s0 (t) + N1 (t); s2 (t) = h2 (t) ∗ s0 (t − ∆τ ) + N2 (t)), (1) where ∆τ is the estimated time delay difference; s1 (t), s2 (t) are two received signals; s0 (t) is the emitted signal; h1 (t), h2 (t) are the impulse responses of media which signals pass; N1 (t), N2 (t) are the noises. In general case, noises of system can be white or colored. The values of si (t), hi (t), s0 (t), Ni (t) can be complex. The time delay estimation methods can be divided in two groups in the frequency and time domain. Methods of estimation in frequency domain consist in the calculation of difference values of phase spectrum, which correspond to the searched signals. Many authors notice the drawbacks of such methods which are 91 connected with influence of the side lobes [2], [5]. Such influence in the frequency domain leads to the shift of estimated value. This effect can be reduced by using the window functions. Other drawback is the problem of measuring values between frequency grid nodes [5]. The time domain delay difference estimation methods can be separated as follows. Methods based on the zero- or threshold- crossings detection. Methods based on calculation of signal model parameters. Methods that use a prior in- formation about phase of signal [6], [7]. Methods based on the cross correlation function analysis [2], and etc. Many authors notice that the phase difference based methods are the most accurate [2], [3]. One of the most effective ways for increasing the accuracy of time delay measurements in the short range radar problem is using the continuous waves with complex modulation. Such signals allow one to provide signal to noise ratio (SNR) higher than pulsed systems with equal power [8]. In radar problems signals with linear frequency modulated continuous waves (FMCW) are often used. One of the main advantages of those one is straight dependence of processed frequency on delay of signals received from each aim. The beat tone obtained by correlation scheme has the following expression [8]: ∆f τ τ 2 ∆f s = A cos(2π[ ]t+2π[f0 τ + +∆θ0 ])+N (t) = A cos(ωt+ϕ))+N (t), (2) Tm 2Tm where A is the amplitude of beat signal; τ is the delay; ∆f is the frequency deviation; f0 is the initial frequency; Tm is the period of modulation; ∆θ0 is the initial phase difference between emitted and received signals; ω is the frequency of the beat tone; ϕ is the initial phase of the beat tone. The beat tone contains information about the time delay of the received signal in its frequency and phase. As a rule, only a frequency is measured [8]. However, it can be shown that the initial phase can also be considered as the information parameter if measured delay difference is less than one period of the signal. Thus, the measured delay τ can be estimated as a function of the frequency and initial phase. It should be noted that in practice the initial phase of beat signal can be considered as a linear function ϕ ≈ 2π[f0 τ +ϕ0 ] if τ  2Tm . The aim of this paper is investigation of possibility of using the phase infor- mation of beat signals for increasing the accuracy of small (less than one period) time delay difference estimation. 2 Problem formulation The cross correlation based method for the time delay estimation is widely used in practice. In one of it implementations, it consists in phase difference estima- tion by using of maximum of the corss-correlation function [5]. The additional advantage of this method is the relatively low computational complexity. It can be shown that in the case of single harmonics the method can be considered as follows [2]: X .qX X  θ̂(t) = arg ρ(0) = arg s2 (t)s∗1 (t) s22 (t) s21 (t) , (3) 92 where ρ(0) is the normalized cross-correlation coefficient; ∗ is conjugation. The advantage of estimator (3) is the relatively low computational complex- ity. However, it is shown that the estimator (3) is not asymptotically effective and has a bias in general case [2].Other drawback of this method is the deteriora- tion in accuracy in areas near 0 and π values of analyzed signals phase difference [1]. The Author of paper [6] has proposed to use a prior information of phase to time relation for parameter estimation of processed signals. The initial phase and frequency of single harmonic signal can be calculated by approximation of full phase as follows: ϕ(t) = ω̂t + θ̂, (4) where ϕ(t) is the approximation line; ω̂ is the estimated frequency; θ̂ is the estimated initial phase. By using the least-square method (LSM ), the initial phase of (4) can be found as [6] (N −1)/2 1 X θ̂ = arg s(n), (5) N n=−(N −1)/2 where N is the sample length; arg is the operator of obtaining phase of a complex number; s(n) is the analyzed sample. In the considered problem s(n) = s2 (n)s∗1 (n), where s2 (n), s1 (n) are complex signals, the delay difference between which is measured. It is noted that solution (5) has the variance which coincides with the Cramer- Rao low bound (CRLB). The Authors of [7] has performed analysis of phase approximation task with assumption that |s(n)| is a random variable. In this case the solution of least square equation has follow expression: N X −1 N X −1 θ̂ = Ψ −1 (−β[ n|s(n)| arg s(n)] + η[ |s(n)| arg s(n)]), (6) n=0 n=0 PN −1 PN −1 PN −1 where Ψ = αη − β 2 ; α = n=0 |s(n)|; β = n=0 n|s(n)|; η = n=0 n2 |s(n)|. It is noted that if |s(n)| =const, (6) is equivalent to (5). The model of real-valued harmonic signals of FMCW radars and correlation scheme of processing (2) can be described as follow: s1 (t) = A1 cos(ω1 t + θ1 ); s2 (t) = A2 cos(ω2 t + θ2 ), (7) where ω1 , ω2 are the frequencies of analyzed beat tones; A1 , A2 , θ1 , θ2 are the corresponding amplitudes and initial phases. The result of multiplication of the signals s1 and s2 (s(n) = s2 (n)s∗1 (n)) without noises N (t) in anaclitic form can be written as p p X ∆f n X s(n) = A0 exp(i2π[ τk + f0 τk ]) = A0 exp(iWτ (n)τk ), (8) Tm fs k=1 k=1 93 where s(n) is the processed signal; ∆ω is the beat frequency difference; ∆θ is the initial phase difference; Σω, Σθ are the frequency and initial phase sum respectively; τk is the delay of the k-th signals; Wτ (n) is the weight coefficient which is given as ∆f n ∆f Wτ (n) = 2π[ + f0 ] = 2π[n + f0 ]. (9) Tm fs N Taking into account the noise term N (t) in case of the single harmonic signal (p = 1) (9) can be given as s = s2 s∗1 = A0 exp iWτ (n)τ + N (n) = A0 [µn + iνn ], (10) where µn = cos [Wτ (n)τ ] + Re[N (n)] and νt = sin [Wτ (n)τ ] + Im[N (n)]. It can be supposed that signal (10) has two random parameters A0 and τ . The Fisher matrix in this case has rank 2. Each diagonal element of the Fisher Matrix has the following expression [9] N −1 N −1 1 X ∂ 2 µn ∂ 2 νn A20 X 2 Jτ τ = [ + ] = W (n), (11) σ 2 n=0 ∂τ 2 ∂τ 2 σ 2 n=0 τ where Jτ τ is the diagonal element of the Fisher matrix. The Cramer-Rao low bound for estimation of parameter τ of signals in from (10) is given by SN R−1 SN R−1 var[τ ]τ = PN −1 ≈ , (12) 2 n=0 Wτ (n) 4π N (f0 + ∆f /2)2 2 where var[τ ]τ is the variance of delay estimation by delay as parameter. The CRLB for estimation of delay by frequency is given by [9] Tm 2 fs 2 12SN R−1 3SN R−1 var[τ ]ω = ( ) ( ) ≈ , (13) ∆f 2π N3 π 2 ∆f 2 N where var[τ ]ω is the variance of delay estimation by frequency as parameter of the signal. In (13), it is supposed that the sample length coincides with one period of modulation (Tm ), and, hence, Tm = N/fs . The CRLB for estimation of the delay by the initial phase is given by [9] 1 2 12SN R−1 3SN R−1 var[τ ]θ = ( ) ≈ 2 2 , (14) 2πf0 N π f0 N where var[τ ]θ is the variance of delay estimation by frequency. From comparing expressions (12)-(14) the following relations can be given (f0 + ∆f /2)2 (f0 + ∆f /2)2 var[τ ]τ ≈ 12 var[τ ]θ ≈ 12 var[τ ]ω . (15) f02 ∆f 2 Relations (15) show that the proposed approach of using delay as parameter gives the minimal CRLB. For instance, if f0 = 2∆f , then expression (15) gives: var[τ ]τ = 18.75var[τ ]θ = 75var[τ ]ω . 94 The results mentioned above show the advantages of using full information of the phase of beat tone (that equivalent to use the delay as parameter) for time delay estimation. The main drawback of this approach is restriction, which is connected with the phase uncertainty if its value crosses ±π bound. Results described above lead to the supposition that it is actual to further investigation of the effective estimator design for the problem of time delay difference measurement of beat signals of FMCW radar systems. 3 The estimator design The beat tone (10) can be transformed in a Fourier sequence with the follow coefficients Cl (n) = exp −iWτ (n)τl , (16) where τl is the value of grid of discreet transform, k = 0, . . . , N − 1 are numbers of sample. The Fourier transform by coefficients (16) is given by N X −1 N X −1 S(l) = S(τl ) = s(n)Cl (n) = s(n)e−iWτ (n)τl , (17) n=0 n=0 where S(l) is the specter by τl grid. It is obvious that the maximum of (17) corresponds to delay. It can be supposed that the analogue of the Nyquist theorem for (17) can be expressed as τN −1 ≤ 2τmax , (18) where τN −1 is the maximum value of the grid τl ; τmax is the maximum allowable delay. Condition (18) is equal to the restriction of the phase in range 2π. Thus, τN −1 ≤ 1/(2f0 ). (19) The width of peak in transform (17) is proportional to 1/∆f . As it shown in [6], the expression (10) can be rewrite as follows p X p X s(n) = s2 (n)s∗1 (n) = A0 exp (i[Wτ (n)τk + νIk (n)]) = |sk (n)|ei arg sk (n) , k=1 k=1 (20) where νIk (n) are the phase noises sk (n). The maximum of transform (17) for signal (20) corresponds to the null of its derivative N −1 p ∂S(τl ) X X |τl =τk = Im[ Wτ (n) |sk (n)|e(i arg sk (n) e−iWτ (n)τl ] = 0, (21) ∂τl n=0 k=0 where Im[] is the imaginary part of complex expression. 95 Follow the approach proposed in [10], the multiplication of exponents in (21) at point τl = τk can be transformed in the Taylor series with restriction to its first term as exp xx→0 ≈ 1 + x. Thus, a solution of (21) can be expressed as N X −1 p X N X −1 p X Wτ (n) |sk (n)| arg sk (n) = Wτ2 (n) |sk (n)|τk , (22) n=0 k=0 n=0 k=0 Analysis of equation (22) for p = 1 gives the following solution PN −1 n=0 Wτ |s(n)| arg s(n) ∆τ̂ = PN −1 2 , (23) n=0 Wτ |s(n)| where ∆τ̂ is the estimated delay difference. In general case, estimator (23) can be expressed as ∆τ̂ = [WτT AWτ ]−1 WτT AΨ, (24) where Wτ = [Wτ (0), . . . , Wτ (N − 1)]T is the coefficient vector; A is the vector of amplitudes of the sample; Ψ is the vector of arguments of the sample; Ψ = [arg(s(0)), . . . , arg(s(N − 1))]T ; A = diag[|s(0)|, . . . , |s(N − 1)|]. In a special case when |s(n)| = const; A =const·diag[1, . . . , 1] the estimator (24) is given as ∆τ̂ = [WτT Wτ ]−1 WτT Ψ. (25) 4 Investigation of the estimator properties The proposed estimator (24) corresponds to expression of the Gauss-Markov theorem, thus this estimator is a linear, unbiased and asymptotically effective [11]. Expression of the variance of the estimator (24) is given as [11] var[τ ] = σ 2 [WτT AWτ ]−1 . (26) In case of (25), the variance coincides with CRLB (12) for Wτ (9). Figure 1 shows the relation of standard deviation (STD) in logarithmic scale to SNR for proposed estimator (prop.τ ) and compared ones (Tretϕ (3); FKϕ (4) and MaxCorr (3)) . Also in Fig. 1 the results are shown for the maximum likelihood (ML) estimator and CRLB (12); (13) and (14). The ML was performed as the procedure of calculation of frequency difference which corresponds to the maximum of spectrums. The spectrum calculation was carried out by the Fast Fourier transform (FFT) with zero padding to 222 sample size. The values that are shown in Fig. 1 were calculated for signals (7) with delays 1 and 1.0001. The beat signals was simulated for the follow configuration of FMCW: initial frequency 100; frequency deviation 50; period of modulation 0.05; sample frequency 10000 (samples size 500 points). All values above are given in relative units. It should be noted that the selected configuration does not interfere with the generality of carried investigation. 96 Fig. 1. Relation of STD in log-scale to SNR for proposed estimator and compared ones The obtained results (see Fig. 1) shows that the proposed estimator has the threshold value of STD in range 6 - 10 dB, which is coincides with the theoretical estimations [12]. The results for other estimators also match with the theoretical suppositions. The values of STD for (24) attain CRLB, which confirms the supposition of it asymptotic effectiveness. Investigation of the influence of a parasitic signal in the analyzed sample on the bias of estimation was performed for proposed estimators (24) and (25) and compared ones. The following model of superposition of fundamental and parasitic beat tones was used: s1 (n) = A1 exp (i[Wτ (n)τ1 ]) + A2 exp (i[Wτ (n)(τ1 + τ21 ]), (27) s2 (n) = A1 exp (i[Wτ (n)(τ1 + ∆τ1 ]) + A2 exp (i[Wτ (n)(τ1 + ∆τ1 + τ21 + ∆τ21 )]), where τ21 = τ2 − τ1 ; ∆τ21 = ∆τ2 − ∆τ1 ; τ1 , τ2 are the delays for the first fundamental and parasitic signals respectively (for s1 (n)); ∆τ1 , ∆τ2 are the delay differences for second and first fundamental and parasitic signals respectively (for s2 (n)). In Figure 2 relations are shown for the relative bias δ in % to the difference of delays (τ21 ) normalized on the τ1 . Obtained values were calculated for the proposed estimators prop. τ1 (24) and prop. τ1a (25) and compared ones for A2 = 0.1A1 and ∆τ12 = 0.00015. The conditions of numerical experiment are similar to those ones that was used for obtaining Fig. 1. The Figure 3 is shown the relations for the most accurate (unbiased) estimators from fig. 2 for δ in range of ±0.2%. The Analysis of the performed results (see Fig. 2) show the advantages of the proposed estimator prop.τa (25). It is also should be noted that ML method has the biggest bias with compared ones. The results of estimators: prop.τ , FK and MaxCorr have similar bias. The real part of the considered model of signals (27) can be transformed in the following manner: s = A1 cos (ϕ1 (t)) + A2 cos (ϕ2 (t)) = A cos (ϕ1 (t) + b), (28) 97 Fig. 2. The relation of bias δ in % to delay difference τ12 normalized on the τ1 for investigated estimators Fig. 3. The relation of bias δ in % to delay difference τ12 normalized on the τ1 for investigated estimators q where A = (A1 + A2 cos ∆ϕ(t))2 + A22 sin2 ∆ϕ(t) and b = arcsin (A2 sin ∆ϕ(t)/A). If A1  A2 and b is sufficiently small than it can be made the assumption that sin b ≈ b then b ≈ A2 sin ∆ϕ(t)/A1 . In this case, the analytic form of the (28) is given as A2 s = A1 exp (i[ϕ1 (t) + sin ∆ϕ(t)]). (29) A1 Using (29) and (27) in (9) leads to the following expression for s(n): A2 s = s2 (n)s1 (n)∗ = |s(n)| exp (i[Wτ (n)∆τ1 + ∆s ]). (30) A1 where ∆s = sin (τ21 + ∆τ21 )Wτ (n) − sin τ21 Wτ (n). After simple transforms ∆s can be expressed as ∆s = 2 sin (∆τ21 Wτ (n)/2) cos (∆τ21 Wτ (n)/2 + τ21 Wτ (n)). (31) 98 It can be supposed that for the most part of applications ∆τ21  1 and ∆τ21 Wτ (n)  1, and, hence, ∆s = ∆τ21 Wτ (n) cos (∆τ21 Wτ (n)/2 + τ21 Wτ (n)). (32) Thus, estimated δ = τ − τ̂ has the following expression PN −1 PN −1 2 A2 n=0 Wτ |s(n)|∆s A2 W (n) cos (∆τ21 Wτ /2 + τ21 Wτ ) δ= PN −1 2 ≈ ∆τ21 n=0 τ PN −1 2 . A1 n=0 Wτ |s(n)| A1 n=0 Wτ (33) The sum in the numerator of (33) can be replaced by the corresponding 3 integral which solution is propositional to the N/(τ21 ∆f ). In the case when τ12  ∆τ12 , the investigation of (33) leads to the follow relation A2 sin τ12 1 δ∼ ∆τ12 2 ∼ 2 ; lim δ ∼ 1/∆f 2 f0 . (34) A1 τ12 τ12 N →∞ The results of the bias relation is the frequency deviation that is normalized by the initial frequency are shown in the Fig. 4. The presented results are ob- tained for δΣ (33), δint the integral replacement (see 34) and calculated biases for prop τa and prop τ . The relation ∆f /f0 changes in range [0-1]. The beat signals configuration is the same as used for 2 obtaining for τ1 = 1; τ21 = 0.1; ∆τ1 = 0.0001 and ∆τ21 = 0.00015. All values are given in relative units. The obtained results confirm ones presented in the Fig. 2 and assumptions (34). Fig. 4. The relation of bias δ in % to frequency deviation normalized on the initial frequency for investigated expressions and proposed estimators. 5 Discussion and conclusion The paper propose a new approach for time delay difference estimation of beat signals, which are formed by heterodyne scheme of FMCW signals processing. The method is based on the idea of using delay as a parameter of received beat tones. The Cramer-Rao low bound analysis for this case shows the advan- tages compared to the traditional parameters: frequency and initial phase. For 99 instance, if f0 = 2∆f , then the variance of considered method is by 19 times smaller then for phase and in 75 times smaller then for frequency as parameters. Based on the proposed method, a new estimator has been designed in class of the smooth parabolic finite difference estimators. It was shown that the con- sidered estimator is a linear, consistent, asymptotically effective, and, also cor- responds to the Gauss-Markov theorem. However, the proposed approach has restrictions, which are connected with cyclic nature of phase. The SNR thresh- old is about 6 dB. Advantages of the developed estimator are shown in comparison with other techniques in such class in the influence of the parasitic signals on the bias of obtained value. The analysis of bias expression shows its dependence on FMCW signals configuration and on time delay difference between the fundamental and parasitic beat tones. The increasing of the last parameter value leads to the bias reduction as square of its value. Results of the carried investigation shows the advantages of the proposed ap- proach, particulary, of estimator prop. τa for the tasks of the time delay difference measurement of beat signals in comparison with the considered traditional tech- niques in the areas of the variance and influence of the parasitic signals on the bias. References 1. Komarov I.V., Smolskiy S.M.: Fundamentals of Short-range FM Radar. ARTECH HOUSE USA. p. 314 (2003) 2. Lio Y.Z., Zhao B.: Phase-shift correlation method fot accurate phase difference estimation in range fider. Application optic. 54(11), 3470–3477 (2015) 3. Atayants B.A., Davydochkin V.M., Ezerskiy V.V., et al.: Precision systems of FMCW short-range radar for industrial applications. ARTECH HOUSE USA p.360 (2014) 4. Kalmykov Al. A., Ronkin M.V.: Study of methods for increasing the accuracy of FM radar systems measurement. CriMiCo2014, 1171–1172 (2014) 5. Bjorklund S.: A survey and comparison of time-delay estimation methods in linear systems. UniTryck: Linkoping, Sweden, p. 169 (2016) 6. Tretter S. A.: Estimating the frequency of a noisy sinusoid by linear regression. IEEE Trans. Inform. Theory. IT-3 1, 832–835 (1985) 7. Fu P., Kam P.Y.: ML Estimation of the Frequency and Phase in Noise. Proceeding of IEEE Globecom 2006, 1–5 (2006) 8. Skolnik M.: Introduction to radar Systems. The McGraw-Hill Book Company, p. 581 (1981) 9. Rife D., Boorstyn R.: Single-tone parameter estimation from discrete-time obser- vations IEEE Transactions on Information Theory. 20(5), 591-598 (1974) 10. Fitz M.P.: Further Results in the Fast Estimation of a Single Frequency IEEE trans. on communications. 42, 862–864 (1995) 11. Luenberger D.: Optimization by vector space methods. New-York: Wiley, p. 340 (1969) 12. Baggenstoss, P. and Kay, S.: On estimating the angle parameters of an exponential signal at high SNR. IEEE Trans. Signal Process. 39, 1203-1205 (1991)