# How to accurately compute the phase lag between two time series with same sampling frequency.

23 views (last 30 days)

Show older comments

Hi everyone,

My data set consists of two time series. Both the series have a consist time lag but that varies (a bit) with tiem as well. So, i am interested to compute that variation in phase lag between two series.

How my data looks:

What i get after cross-correlation approcah:

Here is the dataset detail:

clear all

clc

s=load("Data_timeseries.csv");

ev_time=datenum(s(:,1),s(:,2),s(:,2),s(:,4),s(:,5),s(:,6)); % time

ser1=s(:, 7); % series 1

ser2=s(:, 8); % series 2

Data is also attached, Thank you!

##### 0 Comments

### Accepted Answer

Chunru
on 7 Oct 2022

Edited: Chunru
on 7 Oct 2022

s=load(websave("Data_timeseries.csv", "https://www.mathworks.com/matlabcentral/answers/uploaded_files/1148175/Data_timeseries.csv"));

ev_time=datenum(s(:,1),s(:,2),s(:,2),s(:,4),s(:,5),s(:,6)); % time

ser1=s(:, 7); % series 1

ser2=s(:, 8); % series 2

t = (0:length(ser1)-1);

ser1env = envelope(ser1, 100, 'peak'); % need envolope detection

yyaxis left

plot(t, ser1, 'r', t, ser1env, 'g');

yyaxis right

plot(t, ser2, 'b')

ns = length(ser1);

lwin = 2000;

overlap = 1000;

i1 = 1;

c = 1;

while true

i2 = i1+lwin-1;

if i2>ns

break

end

% remove mean before xcorr

[r,lags] = xcorr(ser1env(i1:i2)-mean(ser1env(i1:i2)), ser2(i1:i2)-mean(ser2(i1:i2)), 'coeff');

[rmax, imax] = max(r);

res(c, :) = [i1, lags(imax)];

i1 = i1 + (lwin-overlap);

c = c + 1;

end

res

% remove mean before xcorr

[r,lags] = xcorr(ser1env-mean(ser1env), ser2-mean(ser2), 'coeff');

figure

plot(lags, r);

grid on

[rmax, imax] = max(r);

tmax = lags(imax);

tmax, rmax

hold on

plot(tmax, rmax, 'ro')

##### 8 Comments

### More Answers (2)

William Rose
on 7 Oct 2022

I will address your quesiton about phase, but first I have some quesitons. Then I adress phase, below.

Your upper plot does not show x or y axis values. When I plot the data, it does not look like your plot.

s=load("Data_timeseries.csv");

t=datenum(s(:,1),s(:,2),s(:,3),s(:,4),s(:,5),s(:,6)); % time (days)

t=(t-t(1)); %remove initial time offset

ser1=s(:, 7); % series 1

ser2=s(:, 8); % series 2

plot(t,ser1,'-g',t,ser2,'-b')

legend('Ser1','Ser2'); xlabel('Time (days)')

Plot of your data, above, does not look like your plot.

Maybe if I remove the mean value from each:

plot(t,ser1-mean(ser1),'-g',t,ser2-mean(ser2),'-b')

legend('Ser1-mean','Ser2-mean'); xlabel('Time (days)')

The plot of your mean-adjusted data, above, looks more like the plot you shared, but still not the same. The ser1 data in your file is a lot noisier than the green trace in your plot. Please comment on this difference.

How did you compute and plot cross correlation? I am asking because the cross corr plot you shared does not look like what I would expect. A cross correlation like what you showed is what happens if you do not first remove the mean from both signals. See plots below.

[rxy,lags]=xcorr(ser1,ser2,'normalized'); %compute cross correlation

subplot(211), plot(lags,rxy,'-b'); %plot cross correlation

xlabel('Lags (samples)'); ylabel('r_{ser1,ser2}') %add labels

dt=(t(end)-t(1))/(length(t)-1); %sampling interval (days)

subplot(212), plot(lags*dt,rxy,'-b'); %plot cross correlation

xlabel('Lag (days)'); ylabel('r_{ser1,ser2}') %add labels

The only difference in the plots above is the x-axis scaling. Neither one matches your plot. That's why I'm curious about how you made your cross corr plot. Both of the plots above are not the true cross correlation, because the mean was not removed. The true cross correlation is below.

[rxy,lags]=xcorr(ser1-mean(ser1),ser2-mean(ser2),'normalized'); %compute cross corr

figure; subplot(211); plot(lags*dt,rxy,'-b') %plot cross corr

xlim([-1 1]); title('Cross Correlation'); grid on %add title, grid

xlabel('Lag (days)'); ylabel('r_{ser1,ser2}') %add labels

subplot(212); plot(lags*dt,rxy,'-b') %plot cross corr

xlim([0 .5]); title('Cross Correlation'); grid on %add title, grid

xlabel('Lag (days)'); ylabel('r_{ser1,ser2}') %add labels

This is more reasonable and what I would expect.

----------

Now, regarding phase between the signals. Do you want one estimate of the phase difference between ser1 and ser2, based on this recording? Or do you want a new phase difference estimate for every cycle of the waves? How do you plan to use this information?

I assume you want a single phase estimate. Otherwise there would have been no reason for you to plot the cross correlation.

The plots of the signals, above, show that each is roughly sinusoidal. The frequency is about 19.5 cycles in ten days, so f~=19.5/10=1.95.

The cross correlation xcorr(ser1,ser2) has a positive peak at time t=+0.24, evident on the cross corr plot above. This means ser1 lags ser2 by 0.24 seconds.

The phase difference, in radians, is

,

or , where td is the lag, in units of time, and f is the frequency, in cycles per unit time.

##### 3 Comments

William Rose
on 7 Oct 2022

William Rose
on 7 Oct 2022

There is a nice reuslt below, it just takes a while to get to it.

To estimate phase lag of ser1 relative to ser2, once per cycle, you will want to compare the zero crossing time differences. Use either positive or negative crossings. Phase lag estimates more frequent than once per cycle do not make sense. The phase difference, in radians, is

where is the zero crossing time difference, in units of time, and f is frequency in cycles per unit time.

s=load("Data_timeseries.csv");

t=datenum(s(:,1),s(:,2),s(:,3),s(:,4),s(:,5),s(:,6)); % time (days)

t=(t-t(1)); %time (days) (initial offset removed)

ser1=s(:, 7)-mean(s(:,7)); % series 1 minus its mean

ser2=s(:, 8)-mean(s(:,8)); % series 2 minus its mean

To estimate the zero crossing times, you will need to remove drift (high pass filter) and noise (low pass filter). Use filters with zero phase for this, because filters with non-zero phase could lead to errors in the phase lag esitmates. Zero phase filters include Matlb's filtfilt() and non-causal FIR filters that are symmetric about the middle.

Visual inspection of the signals shows 19.5 cycles in 10 days, so the dominant frequency is 1.95 cycles/day. We can do low pass and high pass filtering simultaneously with a bandpass filter. Therefore we make a bandpass with cutoff frequenencies 1/day and 4/day.

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

Plot the signls pre and post filtering to make sure the filtering process is working as desired:

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

The plots above are reassuring. The filtered signals look like what we expect and desire.

Find the indices of the upward zero crossings of the filtered signals with an ingenious function from @Star Strider

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

Plot the results:

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

Compute the differences between the zero crossing times, and convert those difference to phases, and plot them versus time.

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

This looks good. The mean phase lag above is a bit less than 180 degrees. The average phase lag for the entire signal, determined in my earlier answer, using the cross correlation, was 168 degrees. So the cycle-by-cycle results are consitent with the overall result.

##### 4 Comments

William Rose
on 8 Oct 2022

### See Also

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!