# 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?

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: <> • The definition of function sincis missing. Oct 9, 2017 at 6:38 • @OlliNiemitalo Added that to the end above. – Peter K. Oct 9, 2017 at 10:53 • Thanks. The signal package's version of sinc is internal to the package and only takes one argument, so the definition is needed. Oct 9, 2017 at 12:25 • @OlliNiemitalo D'oh! I must have defined it several questions ago and it's still in my workspace..... :-) – Peter K. Oct 9, 2017 at 12:34 • Hmm I played with it some more and I think that sinc function should be something like sin(omega*pi)/(omega*pi) and get rid of N completely. Oct 9, 2017 at 13:39 # Lagrange parabolic estimator The standard Lagrange polynomial parabolic interpolation peak finding formula from Peter's answer, $$p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}$$ has bias as function of the true delay$d$if the cross-correlation peak is that of a critically sampled sinc. If the sampling frequency is increased, the sample values are increasingly dominated by the Taylor parabola of the cross-correlation centered at the peak, and Lagrange interpolation becomes more accurate. The Lagrange parabola adapts well to different peak shapes. # Lagrange parabolic estimator with bias correction The bias as function of true delay$d$in the peak time$p$given by the Lagrange parabola can be corrected to zero bias by a post processing step: $$d = \frac{\sqrt{32p^2 + 1} - 1}{8p}\tag{1}$$ There are no guarantees that the correction works for other peak shapes than critically sampled sinc. # Zero bias parabolic estimator Here is presented a noise-robust estimator of the peak time that is zero-bias as function of$d$if the cross-correlation is a noiseless delayed critically sampled sinc. Do not assume good performance if the sinc is not critically sampled or if you have another peak shape. In the limit of zero white additive noise and with a uniform distribution of the fractional part of$d$the new estimate has a root-mean-square error of about$0.84416936$times the root-mean-square error in the samples of a unit amplitude sinc, according to an error propagation analysis using partial derivatives of the sample errors. This figure is better than that of any other estimator of same functional form (a weighted sum of three successive sample values divided by another weighted sum of those) that has zero bias as function of$d.$The delay estimate is: $$p = \frac{\gamma - \alpha}{\alpha + 2 \beta + \gamma},\tag{2}$$ where$[\alpha,\,\beta,\,\gamma]$are the sample values just before, at, and just after the largest sample, respectively. The property of the estimator being zero bias as function of$d$can be expressed as: $$\begin{array}{c}-2 < d < 2\\ [\alpha, \beta, \gamma] = \left[\operatorname{sinc}(-d-1),\, \operatorname{sinc}(-d),\, \operatorname{sinc}(-d+1)\right]\quad\Rightarrow\quad p = d\end{array}\tag{3}$$ Amplitude scaling of the sinc doesn't affect the zero bias property. The estimator can still be interpreted as finding the zero of the derivative of a parabola such as: $$\frac{1}{2}(α + 2β + γ)p^2 + (α - γ)p,\tag{4}$$ which is really no good for interpolation of the amplitude. There may be another parabola that has the same zero of the derivative but can be used for approximate amplitude interpolation. With a critically sampled delayed sinc as input, white noise sensitivity of the estimator depends on$d:$$$\lim_{\operatorname{RMSE}_{[\alpha,\beta,\gamma]}\to0}\frac{\operatorname{RMSE}_{p}}{\operatorname{RMSE}_{[\alpha,\beta,\gamma]}} \approx \begin{cases}\frac{\sqrt{2}}{2}&\text{if }d=0,\\ 2\sqrt{2}&\text{if }|d|=1,\\ \sqrt{\frac{\pi^2d^2(d + 1)^2(d - 1)^2(3d^2 + 1)}{2\sin^2(\pi d)}}&\text{otherwise},\end{cases}\tag{5}$$ where$\operatorname{RMSE}$stands for root mean square error. Figure 1. White noise sensitivity as function of$d$(plot of Eq. 5). # Comparison A variation of the R script in Peter's answer randomizes also delay lengths in addition to the test signal samples, and produces these errors in the delay estimate for the different estimators: Figure 2. Errors in the delay time estimate for 500 simulation runs, with standard deviation in parenthesis. Turquoise: Lagrange parabola (0.08), red: zero bias parabolic estimator (0.016), black: Lagrange parabola with bias correction (0.004). In this case the lagrange parabolic estimator with bias correction would be the best choice to take, because it gives the least severe errors. The R script that created Fig. 2: library(signal) sinc <- function(omega) { #Taylor approximation near omega=0 return((abs(omega) < 10^-10)*(- (pi*omega)^6/5040 + (pi*omega)^4/120 - (pi*omega)^2/6 + 1) + (abs(omega) >= 10^-10)*sin(pi*omega)/(pi*omega) ) } fractional_delay_filter <- function(delay, length) { eps <- 10^-10 t <- seq(0,length-1) t[t == 0 ] <- eps filter <- sinc(t-delay) return(filter) } # Lagrange parabolic peak_location1 <- function(ykm1,yk,ykp1) { return(0.5*(ykm1 - ykp1)/(ykm1 - 2*yk + ykp1)) } # Robust zero-bias parabolic peak_location2 <- function(ykm1,yk,ykp1) { return((ykp1 - ykm1)/(ykm1 + 2*yk + ykp1)) } # Lagrange parabolic with correction peak_location3 <- function(ykm1,yk,ykp1) { p <- 0.5*(ykm1 - ykp1)/(ykm1 - 2*yk + ykp1) return((sqrt(32*p^2 + 1) - 1)/(8*p)) } estimate_delay <- function(x,xd,locator) { cc <- ccf(xd,x,plot=F) idx <- which.max(cc$acf)
return(locator(cc$acf[idx-1], cc$acf[idx], cc$acf[idx+1]) + cc$lag[idx])
}

Nres <- 500

true_delay_frac <- matrix(, nrow = Nres, ncol = 1)
errors <- matrix(, nrow = Nres, ncol = 3)
for (n in seq(1,Nres))
{
true_delay <- runif(1, 10, 11)
d5p1 <- fractional_delay_filter(true_delay,32)

T <- 1000

x <- rnorm(T,1)
xd <- filter(d5p1, 1, x)

true_delay_frac[n, 1] <- true_delay - 10;
errors[n, 1] <- estimate_delay(x,xd[1:T],peak_location1)-true_delay
errors[n, 2] <- estimate_delay(x,xd[1:T],peak_location2)-true_delay
errors[n, 3] <- estimate_delay(x,xd[1:T],peak_location3)-true_delay
}

plot(errors[,1],col='turquoise', ylab="", ylim=c(-0.2,0.2))
par(new=T)
plot(errors[,2],col='red', ylab="", ylim=c(-0.2,0.2))
par(new=T)
plot(errors[,3],col='black', ylab="Error", ylim=c(-0.2,0.2), xlab="Index")
par(new=F)
title('Delay estimation error')

sqrt(var(errors[,1]))
sqrt(var(errors[,2]))
sqrt(var(errors[,3]))


You can fit a curve to the points around the peak of the crosscorrelation ontained by xcor and find the peak of the fitted curve. Ideally you know the cross correlation function of your signals and you fit that function. For practical purposes a parabola also would do. As a rule of thumb for this approach to work properly the bandwidth of your signal should be smaller than sampling rate/4. The bandwidth limitation is true also for fractional delay filters.

Another approach is to look at the signals in the frequncy domain. You can search the literature for that.

You are looking for the cross-correlation peak under the assumption that your sampled signals are low-pass signals. That means that the actual signals can be reconstructed by convolving with an appropriate sinc function.

So you are looking for something like $(a*\textrm{sinc})*(b*\textrm{sinc})(-)$ which simplifies to $a*b(-)*\textrm{sinc}$ since convolving with multiple sinc functions does not change anything after the first convolution.

So basically you can just do your cross correlation in your sampled domain, then convolve the result with a finer-grained sinc function and search for your maximum. For example, take the central part of your cross correlation, use a window function to reduce artifacts, insert something like 7 zeros between every pair of samples, do an FFT, zero out everything but the central 12.5% of coefficients, transform back. Instead of the 7 zeros between samples you will now have more detailed content suitable for looking for a peak.