Unexpected phase delay via cross-correlation ? Method issue or bad data?

Hi everyone,

I am using the cross-correlation method to measure the phase delay between two-time series. However, sometimes I get unexpected results, e.g. sudden change of phase delay from positive to negative. I have made multiple attempts to fix this issue, but it still persisted. May someone here suggest to me how can I tackle this?

Here is my script:

D=readmatrix('data.csv'); 
tim=datenum(D(:,1), D(:,2), D(:,3), D(:,4), D(:,5), D(:,6));
s1=[tim D(:,7)];
s2=[tim D(:,8)];
A=tim(1)
pre=[]

for i=0:1:8;

    t1_f= addtodate(A, i, 'day')
    t2_f= addtodate(A, i+2, 'day')
    sf1=s1(s1(:,1)>= t1_f & s1(:,1)<= t2_f,:);
    sf2=s2(s2(:,1)>= t1_f & s2(:,1)<= t2_f,:);
    if length(sf1) ~= length(sf2)
        continue;
    end
    S=[sf1 sf2];
    tfa=S(:, 1);
    t=(tfa-tfa(1));  %time (days) (initial offset removed)

    % Pre-processing to refien the time series --------%

    y=S(:, 2)-mean(S(:,2)); % series 1 minus its mean
    x=S(:, 4)-mean(S(:,4)); % series 2 minus its mean
    x=detrend(x);
    y=detrend(y);

    N=length(t);
    dt=(t(end)-t(1))/(N-1); %sampling interval (days)
    fs=1/dt;  fNyq=fs/2;    %sampling frequency, Nyquist frequency (cycles/day)
    f1=1; f2=4;             %bandpass corner frequencies (cycles/day)
    [z,p,k]=butter(4,[f1,f2]/fNyq,'bandpass');  %4th order Butterworth bandpass filter coefficients
    [sos,g]=zp2sos(z,p,k); %sos representation of the filter
    xde=filtfilt(sos,g,x);  %apply zero phase filter to ser1
    yde=filtfilt(sos,g,y);  %apply zero phase filter to ser2


    %--- Time series normalization -----%
    x_max=max(xde)
    x_min=abs(min(xde))
    if x_max > x_min 
        xxx=xde/x_max
    else 
        xxx=xde/x_min
    end  
    y_max=max(yde)
    y_min=abs(min(yde))
    if y_max > y_min 
        yyy=yde/y_max
    else 
        yyy=yde/y_min
    end 

    zci = @(v) find(diff(sign(v))>0);  %define function: returns indices of +ZCs
    izc1=zci(xde);                      %find indices of + zero crossings of y1
    izc2=zci(yde);                      %find indices of + zero crossings of y2

    ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
    ZT1 = ZeroX(t(izc1),xde(izc1),t(izc1+1),xde(izc1+1));
    ZT2 = ZeroX(t(izc2),yde(izc2),t(izc2+1),yde(izc2+1));


    fdom=(length(ZT1)-1)/(ZT1(end)-ZT1(1));

    [rxy,lags]=xcorr(xde,yde, 1440, 'normalized');     %estimate cross correlation
    dt=(tfa(end)-tfa(1))/(length(tfa)-1);             %sampling interval (days)
    tlag=lags*dt;
    [~,idx] = max(rxy) ;        %index of max of cross correlation
    tmaxcc = tlag(idx) ;        %time of max of cross correlation
    ph=2*pi*fdom*tmaxcc;  %average phase difference (radians)
    ph_d = rad2deg(ph);  %phase difference (degrees)
    res=[kk i ph_d max(rxy) fdom]
    pre=[pre,res];

end 

Here is what i get for different windows:

phase_pre = reshape(pre,[],length(pre)/5);
addd=phase_pre';
addd =

    1.0000         0  160.7273    0.8829    1.9134
    1.0000    1.0000  166.3427    0.9135    1.8956
    1.0000    2.0000 -191.8560    0.8727    1.9830
    1.0000    3.0000 -205.8676    0.8109    1.9987
    1.0000    4.0000  159.0571    0.8395    1.9221
    1.0000    5.0000 -201.9391    0.7768    1.9141
    1.0000    6.0000 -191.4146    0.7167    1.8230
    1.0000    7.0000  164.0249    0.7882    1.7975
    1.0000    8.0000 -182.1604    0.8674    1.9327

Example data is also attached.

[Mathieu NOE]

hello

you can get 360° jumps so why not take the complementary value , for example if the value is negative then take 360 + value

[Andi]

@Mathieu NOE, We can but the key point is that i alredy restrict the matching point to 360 then why we have negative values.

[rxy,lags]=xcorr(xde,yde, 360, 'normalized'); 

the output is like this:

         0  160.7273    0.8829    1.9134
    1.0000  166.3427    0.9135    1.8956
    2.0000 -178.4707    0.8505    1.9830
    3.0000  166.8927    0.7562    1.9987
    4.0000  159.0571    0.8395    1.9221
    5.0000  167.0065    0.7534    1.9141
    6.0000  164.0697    0.6749    1.8230
    7.0000  161.7780    0.7874    1.7975
    8.0000 -173.9463    0.8593    1.9327

[Mathieu NOE]

hello @Andi

doc says :

r = xcorr(___,maxlag) limits the lag range from -maxlag to maxlag for either of the previous syntaxes.

so your code simply specify a maximum lag of 360 samples either in one direction or the other

for me this can lead to either positive or negative phase values

again, if phase is neg take the 360 complement.

Answer 1

[Bjorn Gustavsson]

The phase-jump is just the branch-cut of angle (or atan2). As a measure of phase-shift between 2 harmonic signals it is equivalent to say that one is lagging behind the other by a phase-shift of 175 degrees or leading by a phase-shift of -185 (which is never(?) seen due to convention), same thing with a phase-shift of 185 and -175 degrees. If you really need continuous-ish phase read up on unwrapping algorithms or look into the different fex-contributions that does phase unwrapping - there are even 2-D variants. This might be a very difficult problem, or an obnoxiously fiddly problem, or a rather benevolent problem - all depending on the variability of the phases. This is also why we have cyclic colormaps.

HTH

[Bjorn Gustavsson]

Did this help solve your problem?