How to plot the time series data after remove the phase lag?

Hi everyone,

My data set consists of two-time series, I have computed the phase lag value for the data set which is 170-degree. Now, I want to subtract this phase lag value and plot the time series, but have no idea how can I do this. May someone help out me here?

Data and script is attached:

r=readmatrix('data.csv');


figure(1)
x0=10;
y0=10;
width=550;
height=150
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(r(:,1),r(:,2),'-b')
yyaxis right
plot(r(:,1),r(:,3),'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

Answer 1

[Mathieu NOE]

hello

I assume this is the result you are looking for

for me 170° is almost a polarity inversion, so I changed the y to -y and used also xcorr to obtain the best time alignment between the two series

r=readmatrix('data.csv');
t = r(:,1);
x = r(:,2);
y = r(:,3);

dt = mean(diff(t));
[c_a, lag_a] = xcorr(x,-y);
[~, i_a] = max(c_a);
t_a = lag_a(i_a)*dt; % lag SER2 vs SER1 (days)
ty = t+t_a;
ind = ty>=0;
ty = ty(ind); % consider only data for t>=0
yy = -y(ind);

figure(1)
x0=10;
y0=10;
width=550;
height=150
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(t,x,'-b')
yyaxis right
plot(ty,yy,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

[Andi]

@Mathieu NOE thank you for help! Well, i just give an example of 170. Phase delay varies between 50 to 180, so we can expect any value for phase shift. Approach proposed by you works well for 180 or 170 but how about when the phase delay is 50 degree.

[Mathieu NOE]

hello again

ok, there is more generic approach with fft / ifft

seems to me 180° would be better than 170° in this example

r=readmatrix('data.csv');
t = r(:,1);
x = r(:,2);
y = r(:,3);

ydft = fft(y);
% apply phasor rotation
angl = 170; % degrees
phasor = exp(j*angl*pi/180);
yy = ifft(ydft.*phasor,'symmetric');

figure(1)
x0=10;
y0=10;
width=550;
height=150
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(t,x,'-b')
yyaxis right
plot(t,yy,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

[Andi]

@Mathieu NOE thanks again for sharing another approach. But, i am a bit confsie about this, because i already used these function while computing the phase lag.

Here, I am sharing the script along with the data. In this particular example, the phase lag is around 165.

Additionlly, may you share what is the j parametere you used in second approch [phasor = exp(j*angl*pi/180);].

clear all
clc


cd_ev=readmatrix('Data_timeseries.csv'); 

t=datenum(cd_ev(:,1),cd_ev(:,2),cd_ev(:,3),cd_ev(:,4),cd_ev(:,5),cd_ev(:,6));

t=(t-t(1));  %time (days) (initial offset removed)
ser1=cd_ev(:, 7)-mean(cd_ev(:,7)); % series 1 minus its mean
ser2=cd_ev(:, 8)-mean(cd_ev(:,8)); % series 2 minus its mean


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
y1=filtfilt(sos,g,ser1);  %apply zero phase filter to ser1
y2=filtfilt(sos,g,ser2);  %apply zero phase filter to ser2

subplot(211), plot(t,ser1,'-g',t,y1,'-k');
xlabel('Time (days)'); legend('Ser1=raw','y1=filtered'); grid on
subplot(212), plot(t,ser2,'-b',t,y2,'-k');
xlabel('Time (days)'); legend('Ser2=raw','y2=filtered'); grid on


zci = @(v) find(diff(sign(v))>0);  %define function: returns indices of +ZCs
izc1=zci(y1);                      %find indices of + zero crossings of y1
izc2=zci(y2);                      %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),y1(izc1),t(izc1+1),y1(izc1+1));
ZT2 = ZeroX(t(izc2),y2(izc2),t(izc2+1),y2(izc2+1));

figure; subplot(211); 
plot(t,y1,'-b',ZT1,zeros(size(ZT1)),'r+'); 
title('y1=filtered ser1, and + zero crossings'); grid on
subplot(212); plot(t,y2,'-b',ZT2,zeros(size(ZT2)),'r+'); 
title('y2=filtered ser2, and + zero crossings'); grid on
xlabel('Time (days)')


tzcd=ZT1-ZT2(1:end-1);      %zero crossing time differeces (days)
fdom=19.5/10;               %dominant frequency (cycles/day)
phz=360*fdom*tzcd;          %phase lag of y1 relative to y2 (degrees)
figure; plot(ZT2(1:end-1),phz,'-rx'); 
xlabel('Time (days)'); ylabel('Lag (degrees)'); title('Phase Lag v. Time')
ylim([0,360]); set(gca,'ytick',0:90:360); grid on

phase=median(phz)

[Mathieu NOE]

hello again

seems I have done the phase correction on the wrong serie , so thanks to your code above I have corrected that

I also changed j to 1i as the new standard matlab definition of imaginary term;

clc
clearvars

r=readmatrix('data.csv');
t = r(:,1);
x = r(:,2);
y = r(:,3);

xdft = fft(x);
% apply phasor rotation
angl = 165; % degrees
phasor = exp(1i*angl*pi/180);
xx = ifft(xdft.*phasor,'symmetric');

figure(1)
x0=10;
y0=10;
width=550;
height=150
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(t,xx,'-b')
yyaxis right
plot(t,y,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

[Andi]

@Mathieu NOE Thanks for suggetsion. I am still a bit confuse about teh phase shift plot.

Here is the orginal time series plot.

If, we consider the lower panel, and assuem a phase shift of 165 degree, its hard to assuem it looks like same as the one you suggested. For example, our data set have 1440 points per day (that is 360 degree), if we assuem the phase shift is 165, it means if we move the time series by a shift of around 660 points it should gives symmetric results. but not exact ruplicate as your plot shows.

Plus, i try to use the

clear all
clc

r=readmatrix('data.csv');


figure(1)

subplot(211); 
x0=10;
y0=10;
width=550;
height=250
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(r(:,1),r(:,2),'-b')
yyaxis right
plot(r(:,1),r(:,3),'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on
s1=r(:,2);
s2=r(:,3);
t=r(:,1);
s3 = circshift(s1,680);

subplot(212); 
x0=10;
y0=10;
width=550;
height=250;
set(gcf,'position',[x0,y0,width,height])
yyaxis left
plot(t,s1,'-b')
yyaxis right
plot(t,s3,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

But, it skip one cycle. I am not sure what's whats wrong here. May you share your though about this issue.

[Mathieu NOE]

IMHO, a time shift method like this one can only be used if the signals are stable in amplitude and repetitive (sine, square waves or any periodic signals that does not have amplitude variations)

here your signals envelops are not flat , so shifting in time will not give you the best visual results

funny also , the same code as yours gives me a slightly different plot :

[Andi]

Yes, my time series data is highly unstable both in amplitude and sampling rate. Well, I think this difference is only because of the phase delay value. Maybe if we both use the same value for phase delay, the result should be the same. I am trying to test all the available techniques for accurate estimation of the phase lag plus the phase lag plotting to verify the estimated results. Would you like to share your thought about the script I shared for my phase delay estimation based on the zero-crossing approach? I am very grateful for your useful insights.

[script is as below]

clear all
clc
cd_ev=readmatrix('Data_timeseries.csv'); 
t=datenum(cd_ev(:,1),cd_ev(:,2),cd_ev(:,3),cd_ev(:,4),cd_ev(:,5),cd_ev(:,6));
t=(t-t(1));  %time (days) (initial offset removed)
ser1=cd_ev(:, 7)-mean(cd_ev(:,7)); % series 1 minus its mean
ser2=cd_ev(:, 8)-mean(cd_ev(:,8)); % series 2 minus its mean
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
y1=filtfilt(sos,g,ser1);  %apply zero phase filter to ser1
y2=filtfilt(sos,g,ser2);  %apply zero phase filter to ser2
subplot(211), plot(t,ser1,'-g',t,y1,'-k');
xlabel('Time (days)'); legend('Ser1=raw','y1=filtered'); grid on
subplot(212), plot(t,ser2,'-b',t,y2,'-k');
xlabel('Time (days)'); legend('Ser2=raw','y2=filtered'); grid on
zci = @(v) find(diff(sign(v))>0);  %define function: returns indices of +ZCs
izc1=zci(y1);                      %find indices of + zero crossings of y1
izc2=zci(y2);                      %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),y1(izc1),t(izc1+1),y1(izc1+1));
ZT2 = ZeroX(t(izc2),y2(izc2),t(izc2+1),y2(izc2+1));
figure; subplot(211); 
plot(t,y1,'-b',ZT1,zeros(size(ZT1)),'r+'); 
title('y1=filtered ser1, and + zero crossings'); grid on
subplot(212); plot(t,y2,'-b',ZT2,zeros(size(ZT2)),'r+'); 
title('y2=filtered ser2, and + zero crossings'); grid on
xlabel('Time (days)')
tzcd=ZT1-ZT2(1:end-1);      %zero crossing time differeces (days)
fdom=19.5/10;               %dominant frequency (cycles/day)
phz=360*fdom*tzcd;          %phase lag of y1 relative to y2 (degrees)
figure; plot(ZT2(1:end-1),phz,'-rx'); 
xlabel('Time (days)'); ylabel('Lag (degrees)'); title('Phase Lag v. Time')
ylim([0,360]); set(gca,'ytick',0:90:360); grid on
phase=median(phz)

[Mathieu NOE]

can y share the file 'Data_timeseries.csv' ?

[Andi]

sure, .. Please see the attachement.

[Andi]

@Mathieu NOE Here i try to test both method rotation one (propsoed by you) and the time shift for another set of data and here is what i get [Data is also attached]:

Rotation method

r=readmatrix('example.csv');
t=r(:,1);
y1=r(:,2);
y2=r(:,3);
xdft = fft(y2);
lgg=187.83;
phasor = exp(1i*lgg*pi/180);
y3 = ifft(xdft.*phasor,'symmetric');

figure(1);
subplot(211); 
x_f=10;
y_f=10;
width=550;
height=250;
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y2,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

subplot(212); 
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y3,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

Result of rotation method is like

Time shift method

r=readmatrix('example.csv');
t=r(:,1);
y1=r(:,2);
y2=r(:,3);
lgg=187.83;
lgg=lgg*4
lgg=round(lgg/2)
y3 = circshift(y2,lgg);
figure(1);
subplot(211); 
x_f=10;
y_f=10;
width=550;
height=250;
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y2,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

subplot(212); 
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y3,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

Result of timeshift method is like

What's your though about these observations.

[Mathieu NOE]

well I think that the first method works better

I could further improve the code so the phase is computed by xcorr method

clear all
clc

r=readmatrix('example.csv');
t=r(:,1);
y1=r(:,2);
y2=r(:,3);

% lag computation 
dt = mean(diff(t));
[c_a, lag_a] = xcorr(y2,y1);
[~, i_a] = max(c_a);
t_a = lag_a(i_a)*dt; % lag SER2 vs SER1 (days)

% period of ser1 & ser2 => computation of phase delta
ZT1 = find_zc(t,y1,0);
ZT2 = find_zc(t,y2,0);
fdom1 = 1/mean(diff(ZT1));
fdom2 = 1/mean(diff(ZT2));
fdom = 0.5*(fdom1+fdom2);
phz=360*fdom*t_a;          %phase lag of y1 relative to y2 (degrees)

xdft = fft(y2);
phasor = exp(1i*phz*pi/180);
y3 = ifft(xdft.*phasor,'symmetric');
figure(1);
subplot(211); 
x_f=10;
y_f=10;
width=550;
height=250;
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y2,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on
subplot(212); 
set(gcf,'position',[x_f,y_f,width,height])
yyaxis left
plot(t,y1,'-b')
yyaxis right
plot(t,y3,'-r');
xlabel('Time (days)'); legend('SER1','SER2'); grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

the results seems quite good, I don't think we can get better results than that !

[Andi]

@Mathieu NOE Thank you very much for your thoughtful suggestions. In the last version of the code, I define limits for zero crossing because most of the time the number of zero-crossings is not equal, and then we need to decide which one we should delete to equalize the number of zero-crossing points. But in your version, there is no option for the such a particular scenario.

[Mathieu NOE]

ok understood

all the best