How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?
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 &lt;- function(delay, length) 
{
  eps &lt;- 10^-10
  t &lt;- seq(0,length-1)
  t[t == 0 ] &lt;- eps
  filter &lt;- sinc(t-delay,length)
  return(filter)
}

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


estimate_delay &lt;- function(x,xd)
{
  cc &lt;- ccf(xd,x)
  idx &lt;- which.max(cc$acf)
  return(peak_location(cc$acf[idx-1], cc$acf[idx], cc$acf[idx+1]) + cc$lag[idx])
}

true_delay &lt;- 10.5
d5p1 &lt;- fractional_delay_filter(true_delay,32)
fz &lt;- freqz(d5p1)
#plot(fz)
gd &lt;- grpdelay(d5p1)
#plot(gd)

T &lt;- 1000
Nres &lt;- 100

results &lt;- rep(0,Nres)
for (n in seq(1,Nres))
{
  x &lt;- rnorm(T,1)
  xd &lt;- filter(d5p1, 1, x)

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

  results[n] &lt;- 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 &lt;- function(omega, N)
{
  eps &lt;- 10^-56
  return( sin(omega*N/2)/(omega/2 + eps) )
}




So changing that sinc function to:

sinc_pjk &lt;- function(x)
{
  return(sin(pi * x) / (pi * x))
}


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

&lt;>