# Find time lag between two time series

I am trying to find the time lag between two time series over t = [0,1000] using MATLAB (not that it matters). The first time series is simply t^2. The second is (t-15)^2 which is, of course, shifted to the right 15 units (e.g., seconds). My approach has been to find the cross correlation (computed using FFT) and then use the maximum of these to determine the appropriate time shift. For some reason, the answer I keep getting is 0, which doesn't make any sense. Can anyone tell me what I'm doing wrong? Thanks!

I've calculated the cross correlation using FFT as shown (all code is in MATLAB; note fft is the Fast Fourier Transform function and ifft is the inverse Fast Fourier Transform function):

t = [1:1000];
a = t.^2;
b = (t-15).^2;

if (length(a) < length(b))
c = [zeros(1,length(b)-1) a zeros(1,length(b)-length(a))];
d = [b zeros(1,length(b)-1)];
else
c = [zeros(1,length(a)-1) a];
d = [b zeros(1,length(a)-length(b)+length(a)-1)];
end

c_fft = fft(c);
d_fft = fft(d);

cross_spectral_density = c_fft.*conj(d_fft);
cross_correlations = ifft(cross_spectral_density);

[max_xcorr shift_idx] = max(abs(cross_correlations));
shift = shift_idx - max(length(a),length(b));

• Please include plots of the inputs and the cross-correlation output Sep 16 '15 at 13:51
• So, I expect you are getting a thing that looks like a triangle in your cross_correlations array. Cross correlation doesn't like DC offsets. Your signal starts at 1 and goes up, so you have an offset.
– JRE
Sep 16 '15 at 14:25
• In addition to what @JRE says about the DC offset, I think your signal is "nonstationary" and so correlations won't tend to do quite what you expect.
– Peter K.
Sep 16 '15 at 14:38

I have no good explanation for why this happens, except that perhaps quadratic signals are nonstationary in that their mean is changing (or unbounded).

If I do something to ensure the mean is bounded (and not changing) such as:

a = sin(t.^2);
b = sin((t-15).^2);


or

a = rem(a,100);
b = rem(b,100);


then I get the right answer: 15 samples.

This is what happens when I use the first:

This is what happens when I use the second:

both modified signals generate the same (correct) estimate of delay.

Scilab code to generate it below. I also modified your zero padding etc. because it seems overly complicated.

Code Only Below

// 25867
t = [1:1000];
a = (t.^2);
b = ((t-15).^2);
a = a - fix(a/200)*200;
b = b - fix(b/ 200)*200;

c = [a zeros(1,length(b) - 1)];
d = [b zeros(1,length(a) - 1)];

mu = mean([mean(c) mean(d)])

c_fft = fft(c - mu);
d_fft = fft(d - mu);

cross_spectral_density = c_fft.*conj(d_fft);
cross_correlations = fftshift(ifft(cross_spectral_density));
xc = xcorr(c- mu,d- mu,length(t)-1);

[max_xcorr shift_idx] = max(abs(cross_correlations));
shift = shift_idx - max(length(a),length(b));

clf;
subplot(211);
plot(a);
plot(b,'r');
subplot(212)
plot(cross_correlations);
plot(xc,'g:')
plot(shift_idx,cross_correlations(shift_idx),'r.')
xtitle('Delay is ' + string(shift))