Unexpected phase delay via cross-correlation ? Method issue or bad data?
5 次查看(过去 30 天)
显示 更早的评论
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 个评论
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Bartlett 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!