How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?

I would like to calculate a delay between two signals in MATLAB. Normally, I would do this with a cross-correlation (xcorr), and then calculate the position of the correlation peak, but in this case I know that the delay is smaller than one sampling period. Therefore, the method I propose will yield either $0$ or $1$ sample as a delay.

One solution I can come up with, is to construct an interpolated version of the two signals, and then calculate the delay. However, then I still rely on the sampling period of this upsampled signals. Another possible solution is to calculate the cross-correlation peak for different blocks, and then calculate the mean delay. However, I am not sure whether this solution will converge to the right delay. A third solution I can come up with, is to apply a fractional delay to one of the two signals by means of a sinc reconstruction, and then calculate for which fractional delay the correlation peak is the largest. This requires a lot of calculations however, and I still rely on the fractional delays I apply (so I end up with the same problem as with the upsampling).

Therefore, I am wondering whether there is a more elegant or precise solution?

Answer

This is basically what @hooman suggests: fit a parabola to the three points near the peak of the sample cross-correlation of the data.

Using the formula for $p$ here: $$ p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma} $$ where $\alpha,\beta,$ and $\gamma$ are the values of the sample cross-correlation just before the peak, at the peak, and just after the peak respectively.

I've assumed that the signal is a white noise sequence and that a simple FIR filter can be used to simulate a fractional delay.

The plot below shows the FIR filter in the top plot, and 100 estimates of the fractional delay using this technique.

---

R Code Below

#35112
library(signal)

fractional_delay_filter <- function(delay, length) 
{
  eps <- 10^-10
  t <- seq(0,length-1)
  t[t == 0 ] <- eps
  filter <- sinc(t-delay,length)
  return(filter)

}

# https://ccrma.stanford.edu/~jos/parshl/Peak_Detection_Steps_3.html
peak_location  <- function(ykm1,yk,ykp1)
{
  return(0.5*(ykm1 - ykp1)/(ykm1 - 2*yk + ykp1))
}


estimate_delay <- function(x,xd)

{
  cc <- ccf(xd,x)
  idx <- which.max(cc$acf)
  return(peak_location(cc$acf[idx-1], cc$acf[idx], cc$acf[idx+1]) + cc$lag[idx])
}

true_delay <- 10.5
d5p1 <- fractional_delay_filter(true_delay,32)
fz <- freqz(d5p1)
#plot(fz)

gd <- grpdelay(d5p1)
#plot(gd)

T <- 1000
Nres <- 100

results <- rep(0,Nres)
for (n in seq(1,Nres))
{
  x <- rnorm(T,1)

  xd <- filter(d5p1, 1, x)

  #plot(x, type='l')
  #lines(xd[1:T],col='red')å

  results[n] <- estimate_delay(x,xd[1:T])
}

par(mfrow=c(2,1))
plot(d5p1, type='l', col='blue')

title('Fractional filter')

plot(results,col='blue')
lines(c(1,Nres), c(true_delay,true_delay), col='green')
title('Estimates and true delay')

---

The sinc function is from the signal library (which you've said privately you have trouble installing in R). It's defined as:

sinc <- function(omega, N)
{

  eps <- 10^-56
  return( sin(omega*N/2)/(omega/2 + eps) )
}

---

So changing that sinc function to:

sinc_pjk <- function(x)
{
  return(sin(pi * x) / (pi * x))
}

and the test delay to 10.1 samples, Olli's version is less biased:

<>