xcorr always gives best aligment of 0 (which by visual inspection is wrong)

i have some datacurves, which im trying to align in time. while there is clearly max correlation that can be seen visually (like in image below), corrcoeff(y1,y2) always gives maximum correlation of length(y) i.e. aligment of 0. what i am doing wrong here and how can I find out the index of maximal similarity of two curves? Visually it can be seen that maximum correlation is at aligment ~30000

arrays y1 and y1 are in attachment.

[Paul]

Hi Sven,

You're more likely to get help if you post your data by saving it to a .mat file and using the paperclip icon on the Insert strip and post the code that shows the unexpected result, preferably by copy/pasting it into your question and using hte code formating on the Code strip.

[Matt J]

It is not clear how corrcoeff would allow you to align the two curves. Perhaps you intended to use xcorr.

[Sven Larsen]

yes, my bad! corrected title....

[Sven Larsen]

added y1y2.mat as attachment.

Answer 1

[Matt J]

Visually it can be seen that maximum correlation is at aligment ~30000

I'm afraid not:

load y1y2;

x=1:numel(y1);

y3=interp1(x,y1,x+30000,'linear',0); %y1 shifted by 30000

corr0=y1*y2' %correlation with no shift

corr0 = 1.4636e+09

cor30000=y3*y2' %correlation when shifted by 30000

cor30000 = 1.2789e+09

plot(x,[y1;y2;y3]'); legend y1 y2 y3;
ylim([70,92])

[Sven Larsen]

yes it is but thats not the point. point is that xcorr doesent give ANY meaningful estimation what would be best fit!

I had to do my own algorithm to achieve that. algorithm which adds same kind of delay like xcorr but actual calculation is:

sum(abs(y1.'-y2.'))

one needs to calculate that in every delay and best fit is lowest value of those values.

Weird that matlab doesnt have this kind of curve aligment....

[Matt J]

yes it is but thats not the point. point is that xcorr doesent give ANY meaningful estimation what would be best fit!

I disagree. Zero-shift is the best fit. My previous plot makes clear that the shift of 30000 that you expect gives a very poor match above x=1.7e5.

I had to do my own algorithm to achieve that. algorithm which adds same kind of delay like xcorr but actual calculation is: sum(abs(y1.'-y2.'))

And you claim this gave you the result you wanted? My implementation below says otherwise. It gives zero-shift (t=0), even when initializing the iterations at t=30000. Perhaps you used some kind of extrapolation rule that you did not tell us about.

load y1y2;

x=1:numel(y1);

t=fminsearch(@(t)  lsqfunc(t,x,y1,y2)  ,30000 )

t = 0

function  fval = lsqfunc(t,x,y1,y2)

  y1= interp1(x,y1,x+t,'linear',0);

  fval=sum(abs(y1.'-y2.'));


end

[Sven Larsen]

unfortunately I dont understand your logic how zero shift can be better fit. here is before/after aligment figures (i also added another pair for example):

it is pretty clear that right figures are aligned correctly while left figures are not. Yet left figures are optimal according to xcorr/alignsignals which is clearly incorrect.

there really should be "minimal area between figures" curve align method or corrcoef align method built in matlab.

[Stephen23]

"unfortunately I dont understand your logic how zero shift can be better fit."

Note this very important remark from Matt J: "Perhaps you used some kind of extrapolation rule that you did not tell us about." So far you have not told us how you are extrapolating the data when you add any lag, e.g. some constant value or assuming cyclical data. We cannot guess this important information.

[Matt J]

there really should be "minimal area between figures" curve align method

The right hand figures only agree with "minimal area" error because you've chopped away the upper part of the signals, excluding them from the comparison. In the left hand figure, the signals have lengths of roughly 2e5, but in the right figure, the lengths are only about 1.8e5.

It's not clear why the right hand part of the signal gets lower priority to you. Regardless, it is simply not a standard kind of comparison, and that is why Matlab doesn't offer an off-the-shelf solution for you.

[Sven Larsen]

yes these datapoints are chopped because these parts of signal cannot be compared into anything; there is no datapoints to compare to.

malab doc alignsignals has this kind of example also. and there alignment is done with rising signal method. IMHO there should also be minimal area difference method.

[Matt J]

IMHO there should also be minimal area difference method.

A minimal area difference criterion will always be a local minimum only. The global minima will always be achieved when the shift is so large as to make the signals completely non-overlap. Furthermore, there will be plenty of cases where these extreme global solutions are also the only local solutions, leading to very strange, as well as non-unique, results.

As an example, consider the two sawtooth pulses below. There are two ways to shift y1 so that your "overlap-only" area criterion is minimized. One way is to shift y1 to the left by 1 so that the two signals overlap only at t=0, y=0. The other is to shift y1 to the right by 1 so that both signals overlap only at t=1,y=1. Not only is the solution non-unique, but since they only overlap at one point, it is of questionable value.

t=linspace(-0.5,2,100)';
msk=abs(t-0.5)<=0.5;

y1=(1-t).*msk;
y2=t.*msk;

plot(t,[y1,y2]); axis equal; xlabel t; ylabel y
legend y1 y2

[Paul]

Aren't the figures on the right hand side in this comment also chopping off the points to the left of zero of the shifted red curve?

@Sven Larsen, would you mind posting your code that yields the shifted red curve?

Answer 2

[Matt J]

You could look for an optimal cyclic shift

load y1y2;

N=numel(y1);

r=ifft( fft(y1).*conj(fft(y2)) ,'symmetric');

[~,t]=max(r);
plot([y1;y2;circshift(y1,-t)]')
legend y1 y2 y1-shift

[Sven Larsen]

Thanks, this seems to work and its fast!

it produces weird artefacts tho, because of circle shift. I wonder if your idea could be adjusted so there is no cyclic looping but instead cycled part is omitted from comparison.

Answer 3

[Matt J]

Here's another possible comparison criterion where the signals are pre-normalized so that their minimum values are zero.

load y1y2;

N=numel(y1);


[r,lags]=xcorr(y1-min(y1),y2-min(y2));

subplot(1,2,1)
[~,t]=max(r);
plot([y1;y2; interp1(   (1:N) , y1,(1:N)+lags(t) ) ]')
legend('y1', 'y2', 'y1-shift','location','south')
axis square

subplot(1,2,2); 
plot(lags,r); axis square
xlim([-4e4,+4e4])
ylabel Correlation; xlabel Lag