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

Question

Andi 2022-10-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1833573-how-to-plot-the-time-series-data-after-remove-the-phase-lag

评论： Mathieu NOE 2022-10-26

采纳的回答： Mathieu NOE

data.csv

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2022-10-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1833573-how-to-plot-the-time-series-data-after-remove-the-phase-lag#answer_1082138

在 MATLAB Online 中打开

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Andi 2022-10-25

编辑：Andi 2022-10-25

在 MATLAB Online 中打开

Data_timeseries.csv

@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)

Andi 2022-10-25

编辑：Andi 2022-10-25

在 MATLAB Online 中打开

@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.

Andi 2022-10-25

在 MATLAB Online 中打开

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)

Andi 2022-10-25

在 MATLAB Online 中打开

example.csv

@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 2022-10-25

在 MATLAB Online 中打开

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 2022-10-25

编辑：Andi 2022-10-26

@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.