0
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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)

corr(a,b)

corr(a-mean(a),b-mean(b))

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)))); 

Corr with IFFT FFT

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).

$\endgroup$
  • $\begingroup$ remember, to do either autocorrelation or cross-correlation using the DFT and iDFT, you need to zero-pad your data to twice the original length to prevent time-aliasing. and then, keep in mind that there is a triangular envelope applied to the correlation result. $\endgroup$ – robert bristow-johnson May 24 '20 at 3:39
  • $\begingroup$ the other thing i would do with this is DC block both signals before correlating them. then you can expect your correlation to have both positive and negative values. and then, if there was a polarity reversal that you hadn't modeled, you can look for the maximum magnitude of the cross-correlation result. $\endgroup$ – robert bristow-johnson May 24 '20 at 3:42
  • $\begingroup$ Many thanks to both of you for your answers! @DanBoschen. Do you think that it would be possible to provide the code for your last plot with implementation of corr = ifft(fft(a)fft*(b)). That would help a lot! Also regarding this plot: why is the argmax 17 (and not 4 or -4 for example) Also as I'm not in the signal processing domain, what is DC (Direct current?)? and OP? $\endgroup$ – ecjb May 24 '20 at 9:04
  • $\begingroup$ @ecjb OP is you the Original Poster. DC is a not so good replacement for the mean value but yes comes from “Direct Current”. I will answer your other question with the code in the response. $\endgroup$ – Dan Boschen May 24 '20 at 13:13
  • $\begingroup$ Many thanks @DanBoschen for your excellent and complete response. I think that made me understand the math behind $xcorr$ function. Although I also have a feeling of what is the circular cross correlation, I still fail to grasp the concept/math behind FFT and IFFT. As I saw on you profile that you teach DSP and Python for the Boston Section IEEE, do you have any particular book / ressources to recommend on the matter? $\endgroup$ – ecjb May 24 '20 at 17:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.