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

8 ビュー (過去 30 日間)
Andi 2022 年 11 月 12 日
コメント済み: Bjorn Gustavsson 2022 年 11 月 29 日
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:
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;
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_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);
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.
##### 3 件のコメント1 件の古いコメントを表示1 件の古いコメントを非表示
Andi 2022 年 11 月 14 日
@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 2022 年 11 月 15 日
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.

サインインしてコメントする。

### 回答 (1 件)

Bjorn Gustavsson 2022 年 11 月 15 日
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
##### 1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Bjorn Gustavsson 2022 年 11 月 29 日
Did this help solve your problem?

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeChebyshev についてさらに検索

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by