Get the lag of a signal with respect to the other with cross correlation

Question

I have 2 vector of length 20 each. The first 15 entries of $a$ are exactly the same as the last 15 entries of $b$. There is random noise in the 5 last entries of $a$ as well as the 5 first entries of $b$. I would like to use cross correlation to detect the lag of 5 between those 2 signals and wrote the following code:

%Matlab Code
clear all; close all;
iteration = 128;
suma = 0;
for it = 1:iteration
  offset = 5;
  range = 100;
  vecLength = 20;
  a = randi(range,1,vecLength);
  b = zeros(1,vecLength);
  b(offset:size(b,2)) = a(1:size(b,2)-offset+1);
  noise = randi(range,1,offset-1);
  b(1:offset-1) = noise;
  c = xcorr(a,b);
  %[c,lags] = xcorr(a,b);
  [argvalue, argmax] = max(c);
  argmax = argmax - length(a);
  suma = suma + argmax;
end
lag = suma/iteration

However, I get an average of around -2.0 for the value of lag. Here are the result of 7 runs:

-1.7344
-2.1875
-2.1719
-2.1641
-1.9531
-2.3750
-2.1172

How can I get the back the value of offset = 5 instead?

Answer 1

The offset is actually 4 samples given the code generated since MATLAB starts indexing at sample 1.

The accumulation from the DC offset (the non-zero average value) creates a larger triangular correlation pattern and together with the random component can sometimes overcome what would otherwise be the strongest correlation at an offset of 4. The solution is to remove any average level prior to doing the correlation.

Given the uniform distribution between 1 and 100, the average will be close to 50 (and random itself so will vary from run to run, but the mean of the random process is 50.5). When both sequences align, the correlation due to the offset alone would on average be $20(\bar x ^2)$, where $\bar x$ is the average value. This result would be close to 50,000 in this case. Since fewer and fewer of the samples are included in the correlation process as the lag increases, this value will linearly decrease to 0 with each shift of the samples. Due to the random element each sample will vary from this triangular shape. Similarly the correlation of the sequence itself when aligned would on average be $15(\bar x^2)$. By removing the mean we eliminate the DC offset effect and allow the repeated sequence to dominate the correlation result when those sequences are aligned.

This demonstrated in the plots below plotting the 20 correlation runs individually of the OP's code to show the variability, using

[c,lags] = xcorr(a,b);
plot(lags,c)

and

[c,lags] = xcorr(a-mean(a),b-mean(b));
plot(lags,c)

The other option derived from the cross-correlation theorem is to compute a circular cross-correlation using the FFT given as:

$$Corr = \text{ifft}(\text{fft}(a)\text{fft}^*(b))$$

Since it is a circular correlation, the correlation from the DC offset is the same throughout the correlation and not strongest at zero lag as in the linear correlation so need not be removed in this case.

Here is the general MATLAB code for providing the circular correlation results with FFTs given two sequences of equal length with $N$ samples each. Since the inverse FFT will provide a result starting at time =0 and is circular in time (such that the last sample is also the result for the first negative sample), the fftshift command is used to center time = 0 in the middle of the plot similar to the plots above.

# lags is different if odd or even length
if (mod(N,2)); lags = -(N-1)/2: (N-1)/2; else; lags = -N/2:N/2-1; end
c = fftshift(ifft(fft(a).*conj(fft(b)))); 

As RBJ points out in the comments, this results in time aliasing, as the correlation is effectively done by wrapping around to the beginning rather than correlating to zeros at the ends of the sequence as in the linear correlation that is performed using xcorr. The sequences can be zero-padded (append 20 zeroes to each sequence) prior to taking the FFT to replicate the linear correlation result, but in this case the mean value would need to be removed for the same reasons given above. Given the correlated samples are random noise I actually prefer the circular correlation personally without zero-padding and removing the mean as it provides a constant noise floor in the result (rather than tapering to zero as you can see in the first two plots above).