![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1167028/image.png)
How to plot the time series data after remove the phase lag?
2 次查看(过去 30 天)
显示 更早的评论
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
采纳的回答
Mathieu NOE
2022-10-24
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1167028/image.png)
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 个评论
Andi
2022-10-24
@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
2022-10-24
hello again
ok, there is more generic approach with fft / ifft
seems to me 180° would be better than 170° in this example
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1167073/image.png)
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
2022-10-25
编辑:Andi
2022-10-25
@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
2022-10-25
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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168333/image.png)
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
2022-10-25
编辑:Andi
2022-10-25
@Mathieu NOE Thanks for suggetsion. I am still a bit confuse about teh phase shift plot.
Here is the orginal time series plot.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168343/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168348/image.png)
Mathieu NOE
2022-10-25
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 :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168523/image.png)
Andi
2022-10-25
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
@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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168633/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168638/image.png)
What's your though about these observations.
Mathieu NOE
2022-10-25
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 !
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1168653/image.png)
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)