Phase difference of two continuous signals.

I have two signals, which has constanly changing frequencies. I want to be able to calculate the phase shift of these two signals at discreet time windows. I have been able to calculate the phase shift of the entire spectrum, however, I am unsure if this method applies to a continously changing frequency. Note, I realized after the fact, that I might need to use an external tigger source to keep the initial frequencies of both signals the same. Will this be an issue, or is it possible to calculate the localised phase shift under these conditions? I am using the Phase Difference Measurement add on to calculate the phase shift in mili seconds.

Below is the code I used:

file_path1 = '.csv';
file_path2 = '.csv';

ref = csvread(file_path1, 16, 0); 
meas = csvread(file_path2, 16, 0);

signal_ref = ref(:,2);
signal_meas = meas(:,2);

% stored variables
dt = 1.6e-10;
fs = 1/dt;
n = length(signal1_new_meas);
fn = fs/(n-1);
ff = fs/2;
T = 0 : dt :(n-1)*dt;
T = T';
Freq = 0: fn: ff;

figure (1), clf;
subplot(211);
plot(T, signal_ref)
xlabel('Time (s)')
ylabel('Voltage (V)')
title('Refrence signal')

subplot(212)
plot(T, signal_meas)
xlabel('Time (s)')
ylabel('Voltage (V)')
title('Measurement signal')

FFT_ref = fft(signal_ref);
FFT_ref = FFT_ref(1:n/2);
FFT_ref = (abs(FFT_ref)).^2;

FFT_meas = fft(signal_meas);
FFT_meas = FFT_meas(1:n/2);
FFT_meas = (abs(FFT_meas)).^2;

figure (2), clf;
subplot(211);
plot(Freq, 10*log10(FFT_ref))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
title('Refrence signal')
set(gca,'xlim',[0 200e3])

subplot(212)
plot(Freq, 10*log10(FFT_meas))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
title('Measurement signal')
set(gca,'xlim',[0 200e3])

% time difference measurement
TimeDiff = phdiffmeasure(signal_meas,signal_ref, fs, 'corr');
TimeDiff = 1000*TimeDiff;
% display the result
disp(['Estimated time lag Y->X (Y w.r.t. X) = ' num2str(TimeDiff) ' ms'])
commandwindow
% plot the signals
figure(1)
plot(T, signal_meas, 'b', 'LineWidth', 1.5)
grid on
hold on
plot(T, signal_ref, 'r', 'LineWidth', 1.5)
xlabel('Time, s')
ylabel('Amplitude, a.u.')
title('Two signals with phase difference')
legend('Measurement signal X', 'Reference signal Y')

The the calculated TimeDiff output is calculated to be -0.54902ms.

Would the calculated phase shift be equal to the entire spectrum in this case? Apologies, the .csv file is rather large (about 1.4GB per .csv file), therefore I didnt include it.

Answer 1

[Mathieu NOE]

hello

maybe you could use this code once you have splitted the large signal into smaller chuncks

x = 0:0.1:30;
y1 = sin(x+0.25).*exp(-x/100); % reference signal 
y2 = 0.5*sin(x+0.45).*exp(-x/100);  % signal shifted by 0.2 rad

figure(1),
plot(x,y1,'r',x,y2,'b');

% find crossing points at y threshold = 0
threshold = 0;
[xc1] = find_zc(x,y1,threshold);
[xc2] = find_zc(x,y2,threshold);

% phase shift  = x distance divided by period and multiplied by 2pi
p1 = mean(diff(xc1)); % mean period of signal 1
p2 = mean(diff(xc2)); % mean period of signal 2

phase_shift = mean(xc1-xc2)/p1*2*pi % is the value used in the data generation 

phase_shift = 0.2000

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Francois]

Thank you so much for the help. I will look at this right now and let you know if it resolved the issue. One more question. I don't know if this is the correct method to split the data for calculation, but for now I use a loop function as described below:

window_size = 1024; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
step_size = window_size * overlap_percentage / 100;

% Calculate number of windows
num_windows = floor((length(signal1_new_meas) - window_size) / step_size) + 1;

% Initialize an array to store results
Reference_signal = zeros(window_size, num_windows);
Measurement_signal = zeros(window_size, num_windows);

% Iterate over the signal with sliding window
for i = 1:num_windows
    % Define start and end indices for the current window
    start_index = (i - 1) * step_size + 1;
    end_index = start_index + window_size - 1;
    
    % Extract the current window
    window_meas = signal_meas(start_index:end_index);
    window_ref = signal_ref(start_index:end_index);
   
    Reference_signal(:, i) = window_ref;
    Measurement_signal(:, i) = window_meas;
    
end

Is this the correct way to split the data into chunks? Also, do I need to apply your answer then to all the generated data chunks. For instance, the code you provided applied to the first data chunk would look like this:

% find crossing points at y threshold = 0
threshold = 0;
[xc1] = find_zc(x,y1,threshold);
[xc2] = find_zc(x,y2,threshold);
% phase shift  = x distance divided by period and multiplied by 2pi
p1 = mean(diff(xc1)); % mean period of signal 1
p2 = mean(diff(xc2)); % mean period of signal 2
phase_shift = mean(xc1-xc2)/p1*2*pi % is the value used in the data generation 

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

Then Zx would be the phase shift for the first data chunk?

Apologies, this is all very new to me and I feel overwhelmed.

Once again thank you for the help.

Regards

[Mathieu NOE]

hello

I made a few corrections regarding number of windows, start and end indexes ...

according to your code , if I reduce the overlap , the number of windows increases - which is obviously incorrect . see the code below.

also I wanted to play with phdiffmeasure that I haven't used so far .

Now it seems to give the correct results as soon as the signal frequency does not change too much within each window, but it's hard to quantify how much frequency change is allowed before results start to diverge from truth

see the two dummy signals , they have constant phase shift (0.2 rad) but variable frequency . The squared x term in the sinus is used to creat a variable frequency.

so if the amplitude factor is very low like 0.0001*x.^2 then the output of phdiffmeasure is correct (lower red plot) : it says 0.2 rad accross the entire time range

if I remove one zero (0.001*x.^2) to speed up the frequency sweep , then we start to see a diverging result (output of phdiffmeasure gives 0.15 rad )

I have tried to play with window size and overlap , but there is no significant impact IMHO.

% signal generation
fs = 20;
dt = 1/fs;
x = (0:dt:500)';
signal_ref = sin(x+0.001*x.^2+0.25).*exp(-x/1000); % reference signal 
signal_meas = sin(x+0.001*x.^2+0.45).*exp(-x/1000);  % signal shifted by 0.2 rad

samples = length(signal_ref);

window_size = 1024/4; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1; % your definition (incorrect) 

shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  
num_windows = fix((samples-window_size)/shift +1);

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    
    % phase difference code (with phdiffmeasure)
    phase_shift_rad(ci) = phdiffmeasure(signal_meas,signal_ref, fs, 'corr');
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)

end

xx = x(time_index); % new x axis 

figure(1),
subplot(3,1,1),plot(x,signal_ref)
subplot(3,1,2),plot(x,signal_meas)
subplot(3,1,3),plot(xx,phase_shift_rad,'r');

[Mathieu NOE]

so this is the same test protocol with "my" solution , either with slow or fast frequency changing signals .

Seems to work better as in both situations I got the correct answer (computed phase delta = 0.2 rad as generated into the signals)

now, I admit it's on clean signals, if we have true (noisy) signals, it may lead to others results.

here the code with "fast" changing frequency and constant angle shift (0.2 rad)

% signal generation
fs = 20;
dt = 1/fs;
x = (0:dt:500)';
signal_ref = sin(x+0.01*x.^2+0.25).*exp(-x/1000); % reference signal 
signal_meas = sin(x+0.01*x.^2+0.45).*exp(-x/1000);  % signal shifted by 0.2 rad

samples = length(signal_ref);

window_size = 1024; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);

% phase difference code
% find crossing points at y threshold = 0
threshold = 0;

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
   
    
    % phase difference code
    % find crossing points at y threshold = 0
    [xc1] = find_zc(x,window_ref,threshold);
    [xc2] = find_zc(x,window_meas,threshold);
    m = min(numel(xc1),numel(xc2));
    xc1 = xc1(1:m);
    xc2 = xc2(1:m);
    
    % phase shift  = x distance divided by period and multiplied by 2pi
    p1 = mean(diff(xc1)); % mean period of signal 1
    p2 = mean(diff(xc2)); % mean period of signal 2
    phase_shift_rad(ci) = mean(xc1-xc2)/p1*2*pi; % is the value used in the data generation 
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)

    
end


xx = x(time_index); % new x axis 

figure(1),
subplot(3,1,1),plot(x,signal_ref)
subplot(3,1,2),plot(x,signal_meas)
subplot(3,1,3),plot(xx,phase_shift_rad,'r');ylim([-1 1]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Francois]

Hello again.

Thank you so much for the help I really appreciate it. I seem to understand the answer a little better now. I implemented the code into my new data set, which I made sure is aligned properly using a trigger. I have appended the new files, since they are a lot smaller than the previous data set. Using your code I have been able to calculate the phase shift. However, I want to know, am I supposed to get the same phase shift value for all the windows? Is this the correct method to use if I intend to measure the phase shift of points along a OFDR (Optical Frequency Domain Reflectrometer) fiber sensor? I have tried the cross correlation method, but since I am new to this stuff I don't know if I am preforming the calculation correctly.

[Francois]

Aplogies, I see now your earlier comment on some corrections.

This is what I have been able to do:

file_path1 = 'C:\Users\fjswa\OneDrive\Desktop\trigger\ref_1.csv';
file_path2 = 'C:\Users\fjswa\OneDrive\Desktop\trigger\meas_1.csv';

% Read the CSV file skipping the first 17 rows
ref = csvread(file_path1, 16, 0); % Assuming data starts from B17 to D, adjust as needed
meas = csvread(file_path2, 16, 0);

signal_ref = ref(:,2);
signal_ref_AI = ref(:,3);

signal_meas = meas(:,2);
signal_meas_AI = meas(:,3);

%Stored variables
dt = 3.2e-6;
fs = 1/dt;
ff = fs/2;
n = length(signal_meas);
fn = fs/n;
Freq = 0: fn: ff;
x = 0:dt:dt*(n-1);

figure(1), clf
subplot(211)
plot(x, signal_ref)
xlabel('time(s)')
ylabel('Voltage (V)')
title('Reference signal')

subplot(212)
plot(x, signal_meas)
xlabel('Time (s)')
ylabel('Voltage (V)')
title('Measurement Signal')

figure(2)
plot(x, signal_meas, x, signal_ref)
xlabel('Time (s)')
ylabel('Volatge (V)')

window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers 

% Calculate number of windows
num_windows = fix((n-window_size)/shift +1);
% phase difference code
% find crossing points at y threshold = 0
threshold = 0;

for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,n);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
   
    
    % phase difference code
    % find crossing points at y threshold = 0
    [xc1] = find_zc(x,window_ref,threshold);
    [xc2] = find_zc(x,window_meas,threshold);
    m = min(numel(xc1),numel(xc2));
    xc1 = xc1(1:m);
    xc2 = xc2(1:m);
    
    % phase shift  = x distance divided by period and multiplied by 2pi
    p1 = mean(diff(xc1)); % mean period of signal 1
    p2 = mean(diff(xc2)); % mean period of signal 2
    phase_shift_rad(ci) = mean(xc1-xc2)/p1*2*pi; % is the value used in the data generation 
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)
    
end
xx = x(time_index); % new x axis
figure(3),
subplot(3,1,1),plot(x,signal_ref)
subplot(3,1,2),plot(x,signal_meas)
subplot(3,1,3),plot(xx,phase_shift_rad,'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

It seems that using your help I am able to produce the phase shift of data chunks for the signals. Is there anything I should be aware of when I processes these signals? Signal processing has a steep learning curve but I have been enjoying learning as much as I can from it. Once again, thank you so much for the help.

Kind regards

[Mathieu NOE]

I have a slight problem with your data , they have both the same sampling rate but there is a constant time shift between both time bases (around 9.4110e-07 s) so we must re laign one data set vs the other so they share the same common time axis before moving further

[Francois]

I see. i suspect its because I didn't trigger the oscilloscope properly. Otherwise, is it acceptable to manually shift the data using code?

[Mathieu NOE]

this is my code , with data time realigned

alsoe the angle is forced into -pi to + pi range

seems the signals are very unstationnary so I wonder about what the result really means

% signals
ref = csvread('ref_1.csv', 16, 0); 
meas = csvread('meas_1.csv', 16, 0);
time_ref = ref(:,1);
signal_ref = ref(:,2);
dt = mean(diff(time_ref));

time_meas = meas(:,1);
signal_meas = meas(:,2);
fs = 1/dt;

% time align both signals 
time = time_meas;
signal_ref = interp1(time_ref,signal_ref,time);
clear time_ref time_meas

samples = length(signal_ref);

window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);

% phase difference code
% find crossing points at y threshold = 0
threshold = 0;

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    window_time = time(start_index:end_index);
   
    
    % phase difference code
    % find crossing points at y threshold = 0
    [xc1] = find_zc(window_time,window_ref,threshold);
    [xc2] = find_zc(window_time,window_meas,threshold);
    m = min(numel(xc1),numel(xc2));
    xc1 = xc1(1:m);
    xc2 = xc2(1:m);
    
    % phase shift  = x distance divided by period and multiplied by 2pi
    p1 = mean(diff(xc1)); % mean period of signal 1
    p2 = mean(diff(xc2)); % mean period of signal 2
    phase_shift_rad(ci) = mean(xc1-xc2)/p1*2*pi; % is the value used in the data generation 
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)

    
end


timex = time(time_index); % new x axis 
phase_shift_rad = wrapToPi(phase_shift_rad);

figure(1),
subplot(3,1,1),plot(time,signal_ref)
xlim([min(time) max(time)]);
subplot(3,1,2),plot(time,signal_meas)
xlim([min(time) max(time)]);
subplot(3,1,3),plot(timex,phase_shift_rad,'r');
xlim([min(time) max(time)]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Mathieu NOE]

funny

we get completely different results with phdiffmeasure, but that was also the case with my clean sine waves

% signals
ref = csvread('ref_1.csv', 16, 0); 
meas = csvread('meas_1.csv', 16, 0);
time_ref = ref(:,1);
signal_ref = ref(:,2);
dt = mean(diff(time_ref));

time_meas = meas(:,1);
signal_meas = meas(:,2);
fs = 1/dt;

time = time_meas;
signal_ref = interp1(time_ref,signal_ref,time);
clear time_ref time_meas

samples = length(signal_ref);

window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    
    % phase difference code (with phdiffmeasure)
    phase_shift_rad(ci) = phdiffmeasure(window_meas,window_ref, fs, 'corr');
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)
    
end

timex = time(time_index); % new x axis 
% phase_shift_rad = wrapToPi(phase_shift_rad);

figure(1),
subplot(3,1,1),plot(time,signal_ref)
xlim([min(time) max(time)]);
subplot(3,1,2),plot(time,signal_meas)
xlim([min(time) max(time)]);
subplot(3,1,3),plot(timex,phase_shift_rad,'r');
xlim([min(time) max(time)]);

[Mathieu NOE]

i am also learning on my side

with phdiffmeasure and 'corr' argument the function returns a time difference (you already knew it)

I used a equivalent code from another project and both results are exactly the same .

% signals
ref = csvread('ref_1.csv', 16, 0); 
meas = csvread('meas_1.csv', 16, 0);
time_ref = ref(:,1);
signal_ref = ref(:,2);
dt = mean(diff(time_ref));

time_meas = meas(:,1);
signal_meas = meas(:,2);
fs = 1/dt;

time = time_meas;
signal_ref = interp1(time_ref,signal_ref,time);
clear time_ref time_meas

samples = length(signal_ref);

window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    
    % phase difference code (with phdiffmeasure)
     time_diff(ci) = phdiffmeasure(window_meas,window_ref, fs, 'corr'); % time difference estimation
    
         % my code with xcorr (time lag)
        [c_a, lag_a] = xcorr(window_ref,window_meas);
        [~, i_a] = max(c_a);
        t_a(ci) = lag_a(i_a)/fs; % lag in seconds
    
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)
    
end

timex = time(time_index); % new x axis 
% phase_shift_rad = wrapToPi(phase_shift_rad);

figure(1),
subplot(3,1,1),plot(time,signal_ref)
xlim([min(time) max(time)]);
subplot(3,1,2),plot(time,signal_meas)
xlim([min(time) max(time)]);
% subplot(3,1,3),plot(timex,time_diff,'r');
subplot(3,1,3),plot(timex,time_diff,'r',timex,t_a,'*k');
xlim([min(time) max(time)]);

but here we have only a time difference, is it what you wanted or a phase difference ? this assumes we know what the signal frequencies are

for a phase output we need to use phdiffmeasure with 'dft' argument

[Mathieu NOE]

Results obtained with phdiffmeasure and 'dft' argument seems at first glance not too bad ...

clc
clearvars


% signals
ref = csvread('ref_1.csv', 16, 0); 
meas = csvread('meas_1.csv', 16, 0);
time_ref = ref(:,1);
signal_ref = ref(:,2);
dt = mean(diff(time_ref));

time_meas = meas(:,1);
signal_meas = meas(:,2);
fs = 1/dt;

time = time_meas;
signal_ref = interp1(time_ref,signal_ref,time);
clear time_ref time_meas

samples = length(signal_ref);

window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  

% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);

% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    
    % phase difference code (with phdiffmeasure)
     phase_shift_rad(ci) = phdiffmeasure(window_meas,window_ref, fs, 'dft'); % % phase difference estimation (rad)
    
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)
    
end

timex = time(time_index); % new x axis 
% phase_shift_rad = wrapToPi(phase_shift_rad);

figure(1),
subplot(3,1,1),plot(time,signal_ref)
xlim([min(time) max(time)]);
subplot(3,1,2),plot(time,signal_meas)
xlim([min(time) max(time)]);
subplot(3,1,3),plot(timex,phase_shift_rad,'r');
xlim([min(time) max(time)]);

[Francois]

Hey I was checking the code on my side now. It doesnt matter if I get the phase shift or the time difference, since they are related. But I suspect your code you provided is correct. I will play around a bit more with the data and the code and let you know how it is looking. Thanks so much for the help so far.

[Mathieu NOE]

ok

[Mathieu NOE]

Judt FYI, I made a slight modification in the code to make the common time vector process a bit more robust

here tested for both initial csv files (ref_1.csv,meas_1.csv) and the 2 last ones (ref trigger start_000_ALL.csv,meas trigger start_000_ALL.csv)

in the last case , I got this result , but as I am at the very beginning of the learning curve regarding OFDR, I can't really interpret the results obtained here - hope you can do better than I !

% Results obtained with   phdiffmeasure and 'dft' argument 

% signals
% ref = csvread('ref_1.csv', 16, 0); 
% meas = csvread('meas_1.csv', 16, 0);
ref = csvread('ref trigger start_000_ALL.csv', 16, 0); 
meas = csvread('meas trigger start_000_ALL.csv', 16, 0);
time_ref = ref(:,1);
signal_ref = ref(:,2);
dt = mean(diff(time_ref));
time_meas = meas(:,1);
signal_meas = meas(:,2);
fs = 1/dt;

% define which data has time delay and interpolate the other data on
% correct time data
time_ref_min = time_ref(1);
time_meas_min = time_meas(1);

if time_ref_min>time_meas_min
    time = time_ref;
    signal_meas = interp1(time_meas,signal_meas,time);
else
    time = time_meas;
    signal_ref = interp1(time_ref,signal_ref,time);
end
clear time_ref time_meas

samples = length(signal_ref);
window_size = 32; % Amount of data point used to seperate data chunks.
overlap_percentage = 50; % Desired overlap percentage
overlap = window_size * overlap_percentage / 100;
shift = window_size-overlap;    % nb of samples between 2 contiguous buffers  
% Calculate number of windows
% num_windows = floor((m - window_size) / overlap) + 1;
num_windows = fix((samples-window_size)/shift +1);
% Iterate over the signal with sliding window
for ci = 1:num_windows
    % Define start and end indices for the current window
    start_index = 1+(ci-1)*shift;
    end_index = min(start_index+ window_size-1,samples);
    
    % Extract the current window
    window_ref = signal_ref(start_index:end_index);
    window_meas = signal_meas(start_index:end_index);
    
    % phase difference code (with phdiffmeasure)
     phase_shift_rad(ci) = phdiffmeasure(window_meas,window_ref, fs, 'dft'); % % phase difference estimation (rad)
    
    time_index(ci) = round((start_index+end_index)/2);  % time index expressed as sample unit (dt = 1 in this simulation)
    
end
timex = time(time_index); % new x axis 
% phase_shift_rad = wrapToPi(phase_shift_rad);
figure(1),
subplot(3,1,1),plot(time,signal_ref)
xlim([min(time) max(time)]);
subplot(3,1,2),plot(time,signal_meas)
xlim([min(time) max(time)]);
subplot(3,1,3),plot(timex,phase_shift_rad,'r');
xlim([min(time) max(time)]);