3
$\begingroup$

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));
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Please include plots of the inputs and the cross-correlation output $\endgroup$ – endolith Sep 16 '15 at 13:51
  • 1
    $\begingroup$ 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. $\endgroup$ – JRE Sep 16 '15 at 14:25
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter K.♦ Sep 16 '15 at 14:38
2
$\begingroup$

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:

enter image description here

This is what happens when I use the second:

enter image description here

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))
Improve this answer
$\endgroup$

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.