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

5 次查看（过去 30 天）
Andi
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:
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.
  3 个评论
显示 1更早的评论
Andi
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
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
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

类别

在 Help Center 和 File Exchange 中查找有关 Chebyshev 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by