Estimate Delay of a Known Signal Delayed by Sub Sample Resolution

Question

Given a known signal $ x \left( t \right) $ and its delayed version $ y \left(t, \tau \right) = x \left( t - \tau \right) $.
Both are sampled by Sampling Frequency $ {F}_{s} $ to generate the signals $ x \left[ n \right] $ and $ y \left[ n \right] $.

Assume they are sampled on a time support much larger than their non zero time (Namely, the whole part of the signal where it is not vanished is sampled).

Given those 2 sampled signals, using Cross Correlation (Matched Filter) the Delay Parameter $ \tau $ can be estimated.

Yet, its estimation is limited to the grid defined by the sampling Frequency which means if $ \tau < {T}_{s} = \frac{1}{{F}_{s}} $ the estimation isn't accurate (Actually, the part which is the modulo of the time interval can't be estimated).

Is there an effective method to estimate the delay in Sub Sampling resolution?

Thank You.

Answer 1

Since the Matched Filter is used the Cross Correlation becomes "Auto Correlation" which is assured to be "Symmetric" relative to its maximum.

Hence a good approximation would be to approximate the area around the sampled peak by a Symmetric Function as a parametric model of the function.

The most trivial and easiest model for such function would be the Parabolic Function which is illustrated in the following.

Here's some scilab code that does the cross-correlation and then finds the inter-sample peak by doing a simple (ish) parabolic interpolation.

The plot shows the absolute relative error between the truth and the estimate.

// 25684
vals= [];
taus = [-0.5:.005:0.5] + 10^-9; // Add a small amount to not have a zero.

for tau = taus

    N = 1000;
    t = 0:N-1;
    omega=0.1902834909;
    phi = 2*%pi*0.3823483924;

    x = sin(omega*t + phi);
    x_tau = sin(omega*(t-tau) + phi);

    xc = xcorr(x_tau,x);

    [mx,ix] = max(xc);

    x_vals = [-1 0 +1];
    y_vals = xc(ix + x_vals);

    //https://stackoverflow.com/a/717791/12570
    den = (x_vals(1) - x_vals(2))*(x_vals(1) - x_vals(3))*(x_vals(2) - x_vals(3));
    A = x_vals(3)*(y_vals(2) - y_vals(1)) + x_vals(2)*(y_vals(1) - y_vals(3)) + x_vals(1)*(y_vals(3) - y_vals(2))
    A = A/den;

    B = x_vals(3)^2*(y_vals(1) - y_vals(2)) + x_vals(2)^2*(y_vals(3) - y_vals(1)) + x_vals(1)^2*(y_vals(2) - y_vals(3));
    B = B/den;

    C = x_vals(2) * x_vals(3) * (x_vals(2) - x_vals(3))* y_vals(1) + x_vals(3) * x_vals(1) * (x_vals(3) - x_vals(1))* y_vals(2) + x_vals(1) * x_vals(2) * (x_vals(1) - x_vals(2))* y_vals(3);
    C = C/den;

    x_v = -B / (2*A);
    y_v = C - B*B / (4*A);

    vals= [vals; tau x_v y_v];

end 

clf
plot(taus,abs((vals(:,2)' - taus)./taus));
xlabel('True value of delay');
title('Error in using parabola peak as subsample delay estimate');

Answer 2

Because a shift in time manifests itself as a phase angle shift in the frequency domain, perhaps comparing the spectral phase of the DFT of x(n) and the spectral phase of the DFT of y(n) would be useful. Just a thought.

Answer 3

As the question as been focused, I edit the answer in a three-fold way:

One first approach consists in (somehow) compensating the unknow delay, by applying the initial algorithm for several delays, and in comparing results. This can be done brute force or with optimization loops. One frequent filter design for such a task is known as the Thiran Allpass Interpolators. A Matlab implementation is given at: Thiran Allpass Interpolation in Matlab. A reference paper is J.-P. Thiran, Recursive digital filters with maximally flat group delay, 1971, IEEE Transactions on Circuits Theory.

A second approach may consist in finding a more invariant transformation. The phase suggested by @Richard Lyons is a very good option, since several Fourier transforms possess nice shift/time invariant properties. The two above methods can be extented with GCC or Generalized Cross Correlation, a link with nice background references (like Optimum Estimation of Time Delay by a Generalized Correlator, IEEE Trans. Acoustics, Speech and Signal Processing, 1979). They can also be merged, for instance with GCC-PHAT Cross-Correlation, including phase information.

Invariance is indeed a neat concept, and may be used as well using other signal properties (or non-properties). If the process under consideration is quite stationary, you can get good results with the methods above, up to an integer number of periods if the signal is periodic.

In the case the data is unstationary, has amplitude variations, gets noisy, other methods, or combinations, can be useful. Wavelets can deal with unstationarity. Some are quite shift-invariant. Some others are complex and provide local phase information.

I would suggest a combination of continuous complex wavelets and correlation, such as the one proposed in Time Delay Estimation Using Wavelet Transform for Pulsed-Wave Ultrasound, Annals of Biomedical Engineering, 1995. Other reference? On Wavelet Denoising and its Applications to Time Delay Estimation, IEEE Transactions on Signal Processing, 1999.

A careful reference analysis or the above mentioned concepts may provide you with other useful insights.

My personal experience is not global. It is mostly based on image processing experience, and a work on local delay adaptation for adaptive pattern subtraction in seismic data, using complex wavelets and very short one-tap (unary) complex phase filters, which proved to work quite well.

Answer 4

You can interpolate samples at fractional positions by using a high quality interpolation function or kernel, such as a wide enough windowed Sinc kernel.