Dynamic cross correlation for time delay of two signals

I have two signals, both recorded at the same sampling frequency, however one has a phase delay. I am currently playing around trying to figure out how to calculate the time delay between these two signals using the cross correlation method. I want to calculate the time delay between these two signals using windows, as I want the dynamic time delay of the signals and not the time delay over the entire spectrum. How do I do this? I am using matlab and for now I have this code, but I am unsure whether or not I am doing the right thing. How would I be able to get the maximum peaks of the correlation of all the columns in my "corr" matrix, and how do I convert this into a time delay?

file_path1 = '\ref trigger new_000_ALL.csv';
file_path2 = '\meas trigger new_000_ALL.csv';

% Read the CSV file skipping the first 17 rows
ref = csvread(file_path1, 16, 0); 
meas = csvread(file_path2, 16, 0);

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

n = length(signal_ref);
dt = 3.2e-6;
fs = 1/dt;

FFT_ref = fft(signal_ref);
FFT_ref = FFT_ref(1:n/2);
FFT_meas = fft(signal_meas);
FFT_meas = FFT_meas(1:n/2);

FFT_n = length(FFT_meas);

window_size = 46; % Amount of data point used to separate 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((FFT_n-window_size)/shift +1);
corr_values = cell(num_windows, 1);
lags_values = cell(num_windows, 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,FFT_n);
    
    % Extract the current window
    window_ref = FFT_ref(start_index:end_index);
    window_meas = FFT_meas(start_index:end_index);

    ifft_ref = ifft(window_ref);
    ifft_meas = ifft(window_meas);
    
   [corr_temp, lags_temp] = xcorr(ifft_ref, ifft_meas);
    
    corr_values{ci} = corr_temp;
    lags_values{ci} = lags_temp;
    
end

% Concatenate all correlation and lags values
corr = cell2mat(corr_values');
lags = cell2mat(lags_values);

lag = lags(1,:);

figure(1)
plot(Time, signal_ref);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Reference signal')

figure(2)
plot(Time, signal_meas)
xlabel('Time (s)')
ylabel('Voltage (V)')
title('Measured signal')

figure(3)
plot(Freq, real(FFT_ref));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Reference signal')

figure(4)
plot(Freq, real(FFT_meas))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Measured signal')

figure(5)
plot(lag, abs(corr(:,20)));
xlabel('Cross correlation shift (GHz)')
ylabel('Amplitude (A.U.')

[Mathieu NOE]

hello again

I am not sure to see the difference with your other post

Phase difference of two continuous signals. - MATLAB Answers - MATLAB Central (mathworks.com)

[Francois]

Hello again, I didnt want to post too much on the other post since there is already so much. Thank you again for helping me. I was playing around with the updated code you helped me with, but I was unable to see much change. I have a lot of questions since I am fairly confused with this. Like for instance, if I want to calculate the time delay using this method, then I have to look at the maximum peaks and their corresponding lag values. If say the max peak is at -1, does this mean that the time delay is equal to the amount of data points in the window multiplied with the sampling interval? I am supposed to do signal processing for an OFDR (Optical Frequency Domain Reflectometer) fiber sensor, but I don't know how to calculate the strain along the fiber sensor. If I use this method then the length change along the fiber is excessively high, which should not be the case (usually the length change is in the micro - to nano - meters). The calculation is supposed to be dynamic, since I want the strain at every point along the sensor. However, preforming this calculation on the entire spectrum I think might be wrong. But because I have never done this kind of work, I don't know if what I am doing is correct. Apologies for making another post.

Regards.

[Mathieu NOE]

ok, I'd like to try something but I need your data files

all the best

[Francois]

Hello again, here are the files for the code above.

[Mathieu NOE]

hello again

I tested these files with the code from the other post (I answered a few inutes ago) . here I see you are using a different approach . I am not yet sure to understand how this works and what would be the differences with the code in the other post .

Still need to dive a bit deeper in OFDR processing (I have downloaded a few papers but it's no simple matter ! this is probably going to take me a while to have a basic understanding)

[Mathieu NOE]

I continue to play with your code, that I modified also a bit in some places.

• improved , robustified "common" time section + interpolation . Now the time vector is taken between the valid limits (start with the max initial time of both records and stops with the min of both; typically we will not pick the first(s) samples of the data that has some time advance, and same at the end , we truncate the record that stops in latest position)
• a bit rewamped the fft computations (nothing special, but that was triggered by the fact I could not find where "Freq" is generated in the code you posted in first place)
• wonder why we would not use all xcorr outputs and why you picked one only ? see my results , and FWHM result : 1.366929294441836

clc
clearvars
close all

% 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 valid common Time vector and interpolate data if needed (on new common Time vector) 
% (define which data has time delay and interpolate the other data on
% correct Time vector)
time_ref_min = time_ref(1);
time_meas_min = time_meas(1);
time_ref_max = time_ref(end);
time_meas_max = time_meas(end);

if time_ref_min>time_meas_min
    ind = time_ref<=time_meas_max;
    Time = time_ref(ind);
    signal_ref = signal_ref(ind);
    signal_meas = interp1(time_meas,signal_meas,Time);
else
    ind = time_meas<=time_ref_max;
    Time = time_meas(ind);
    signal_meas = signal_meas(ind);
    signal_ref = interp1(time_ref,signal_ref,Time);
end
clear time_ref time_meas

% FFT process
N = length(signal_ref);
Fs = 1/dt;
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
Freq=k/T;    %create the frequency range

FFT_ref = fft(signal_ref); % make up for the lack of 1/N in Matlab FFT
FFT_meas = fft(signal_meas);% make up for the lack of 1/N in Matlab FFT

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2); 
Freq = Freq(1:cutOff);
FFT_ref = FFT_ref(1:cutOff);
FFT_meas = FFT_meas(1:cutOff);


FFT_n = length(FFT_meas);
window_size = 100; % Amount of data point used to separate 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((FFT_n-window_size)/shift +1);
corr_values = cell(num_windows, 1);
lags_values = cell(num_windows, 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,FFT_n);
    
    % Extract the current window
    window_ref = FFT_ref(start_index:end_index);
    window_meas = FFT_meas(start_index:end_index);
    ifft_ref = ifft(window_ref);
    ifft_meas = ifft(window_meas);
    
   [corr_temp, lags_temp] = xcorr(ifft_ref, ifft_meas);
    
    corr_values{ci} = corr_temp;
    lags_values{ci} = lags_temp;
    
end
% Concatenate all correlation and lags values
corr = cell2mat(corr_values');
lags = cell2mat(lags_values);
lag = lags(1,:);

figure(1)
plot(Time, signal_ref,Time, signal_meas);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Time signals')
legend('Reference signal','Measured signal')


figure(2)
plot(Freq, abs(FFT_ref));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Reference signal')

figure(3)
plot(Freq, abs(FFT_meas))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Measured signal')

figure(4)
% plot(lag, abs(corr(:,20)));
m_corr = mean(abs(corr),2);
plot(lag, m_corr);
hold on
xlabel('Cross correlation shift (GHz)')
ylabel('Amplitude (A.U.')
xlim([-20 20])

% FWHM computation (dx of points at 50% of the peak value)
[ZxP,ZxN] = find_zc(lag, m_corr,max(m_corr)/2);
FWHM = ZxN - ZxP;
plot([ZxP  ZxN], ones(2,1)*max(m_corr)/2,'r');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% 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 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
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 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Francois]

Hello, I have made some changes to the code. For every peak in the 'corr' matrix the FWHM is calculated using a loop and interpolation. I am pretty sure this interpolation is done incorrectly. But, when I tried to use your function 'find_zc' the values for the ZxN and ZxP for some columns were not matching. I know there is probably a better way to find the FWHM the way I did it but for my one data set it was working. However, for my other data set it does not want to work. Some problem with the line:

newLag = linspace(leftIndex, rightIndex, Num_points); % New lag values for interpolation

Apparently it is not a scalar value. I am working on trying to fix it, but this is how far I have gotten for my origninal set of data in my first post. I see that I didn't include the line of how I calculate the variable 'Freq' in the start of this post. Apologies. I only showed one 'corr' peak in this post as an example but I do want to use all of them. I should have stated that from the beginning to avoid confusion.

clc
clearvars
close all

% Results obtained with   phdiffmeasure and 'dft' argument 
file_path1 = 'ref_1.csv';
file_path2 = '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);
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 valid common Time vector and interpolate data if needed (on new common Time vector) 
% (define which data has time delay and interpolate the other data on
% correct Time vector)
time_ref_min = time_ref(1);
time_meas_min = time_meas(1);
time_ref_max = time_ref(end);
time_meas_max = time_meas(end);

if time_ref_min>time_meas_min
    ind = time_ref<=time_meas_max;
    Time = time_ref(ind);
    signal_ref = signal_ref(ind);
    signal_meas = interp1(time_meas,signal_meas,Time);
else
    ind = time_meas<=time_ref_max;
    Time = time_meas(ind);
    signal_meas = signal_meas(ind);
    signal_ref = interp1(time_ref,signal_ref,Time);
end
clear time_ref time_meas

% FFT process
N = length(signal_ref);
Fs = 1/dt;
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
Freq=k/T;    %create the frequency range

FFT_ref = fft(signal_ref); % make up for the lack of 1/N in Matlab FFT
FFT_meas = fft(signal_meas);% make up for the lack of 1/N in Matlab FFT

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2); 
Freq = Freq(1:cutOff);
FFT_ref = FFT_ref(1:cutOff);
FFT_meas = FFT_meas(1:cutOff);

FFT_n = length(FFT_meas);
window_size = 16; % Amount of data point used to separate 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((FFT_n-window_size)/shift +1);
corr_values = cell(num_windows, 1);
lags_values = cell(num_windows, 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,FFT_n);
    
    % Extract the current window
    window_ref = FFT_ref(start_index:end_index);
    window_meas = FFT_meas(start_index:end_index);
    ifft_ref = ifft(window_ref);
    ifft_meas = ifft(window_meas);
    
   [corr_temp, lags_temp] = xcorr(ifft_ref, ifft_meas);
    
    corr_values{ci} = corr_temp;
    lags_values{ci} = lags_temp;
    
end
% Concatenate all correlation and lags values
corr = cell2mat(corr_values');
lags = cell2mat(lags_values);
lag = lags(1,:);

figure(1)
plot(Time, signal_ref,Time, signal_meas);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Time signals')
legend('Reference signal','Measured signal')

figure(2)
plot(Freq, abs(FFT_ref));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Reference signal')

figure(3)
plot(Freq, abs(FFT_meas))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Measured signal')

% Calculate the absolute values of 'corr_'
corr_abs = abs(corr);

% Initialize array to store FWHM values
fwhm_values = zeros(1, size(corr_abs, 2));

% Loop over each column in the matrix
for i = 1:size(corr_abs, 2)
    % Extract column data
    y = corr_abs(:, i);
    
    % Find maximum value and its index
    [maxValue, maxIndex] = max(y);
    
    % Determine half maximum value
    halfMaxValue = maxValue / 2;
    
    % Find indices where signal crosses half maximum on both sides
    leftIndex = find(y(1:maxIndex) <= halfMaxValue, 1, 'last');
    rightIndex = find(y(maxIndex:end) <= halfMaxValue, 1) + maxIndex - 1;

    % Interpolation
    N = 64; % Number of points between the leftIndex and rightIndex
    newLag = linspace(leftIndex, rightIndex, N); % New lag values for interpolation

    new_corr = interp1(1:numel(y), y, newLag, 'linear');
    
    % Find the FWHM from the interpolated data
    [~, maxIndex] = max(new_corr);
    leftIndex = find(new_corr(1:maxIndex) <= halfMaxValue, 1, 'last');
    rightIndex = find(new_corr(maxIndex:end) <= halfMaxValue, 1) + maxIndex - 1;
    
    % Calculate FWHM
    fwhm_values(i) = newLag(rightIndex) - newLag(leftIndex);
end

Toa = fwhm_values * 1e9; %convert FWHM values to GHz 
Time_delay = 1./Toa; %Calculate the time delay
L_i = (Time_delay * 3e8)/(2*1.5); %Calculate the length change of each of the windows

F_L = linspace(0, Fs/2,num_windows); %frequency at every window
L_n = 147059/310; % Calculating frequency amount for every meter (310m fiber sensor). (Quick note, it would be preffered if I just used the frequency bin associated with the max value in the Frequency domain of the signal, currently working to improve this)
Theory_L = F_L./L_n; %Calculate the theoretical distance traveled to each window.

Length_change = Theory_L - L_i; %Calculate the change in the length of the fiber at each window
strain = Length_change./Theory_L; %Calculating the strain over the entire length of the sensor

figure(4)
plot(Theory_L, strain)
xlabel('Distance (m)')
ylabel('Strain (a.u.)')
title('Strain over the length of the fiber')

This graph looks like I am on the right track. However, like I said earlier my loop function to calculate the FWHM and interpolation is most likely wrong. Resulting in a graph like this for the data set named : ref trigger start_000_ALL.csv and meas trigger start_000_ALL.csv.

I see that for this data set the 'fwhm_values' have some 0 values. Is there any reason why this happens?

Thanks again for the help.

[Mathieu NOE]

hmm , ok, I will look at your new code a bit latter (I am a bit busy today)

but just to prove that the FWHM computaion with find_zc function works, here my last code a bit expanded to plot each individual corr FMWH values

see, it works on all data , just I needed to restraint the ZxP,ZxN values to the ones closest to the peak.

that's why you had the issues (should be OK now).

clc
clearvars
close all

% 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 valid common Time vector and interpolate data if needed (on new common Time vector) 
% (define which data has time delay and interpolate the other data on
% correct Time vector)
time_ref_min = time_ref(1);
time_meas_min = time_meas(1);
time_ref_max = time_ref(end);
time_meas_max = time_meas(end);

if time_ref_min>time_meas_min
    ind = time_ref<=time_meas_max;
    Time = time_ref(ind);
    signal_ref = signal_ref(ind);
    signal_meas = interp1(time_meas,signal_meas,Time);
else
    ind = time_meas<=time_ref_max;
    Time = time_meas(ind);
    signal_meas = signal_meas(ind);
    signal_ref = interp1(time_ref,signal_ref,Time);
end
clear time_ref time_meas

% FFT process
N = length(signal_ref);
Fs = 1/dt;
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
Freq=k/T;    %create the frequency range

FFT_ref = fft(signal_ref); % make up for the lack of 1/N in Matlab FFT
FFT_meas = fft(signal_meas);% make up for the lack of 1/N in Matlab FFT

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2); 
Freq = Freq(1:cutOff);
FFT_ref = FFT_ref(1:cutOff);
FFT_meas = FFT_meas(1:cutOff);


FFT_n = length(FFT_meas);
window_size = 100; % Amount of data point used to separate 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((FFT_n-window_size)/shift +1);
corr_values = cell(num_windows, 1);
lags_values = cell(num_windows, 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,FFT_n);
    
    % Extract the current window
    window_ref = FFT_ref(start_index:end_index);
    window_meas = FFT_meas(start_index:end_index);
    ifft_ref = ifft(window_ref);
    ifft_meas = ifft(window_meas);
    
   [corr_temp, lags_temp] = xcorr(ifft_ref, ifft_meas);
    
    corr_values{ci} = corr_temp;
    lags_values{ci} = lags_temp;
    
end
% Concatenate all correlation and lags values
corr = cell2mat(corr_values');
lags = cell2mat(lags_values);
lag = lags(1,:);

figure(1)
plot(Time, signal_ref,Time, signal_meas);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Time signals')
legend('Reference signal','Measured signal')


figure(2)
plot(Freq, abs(FFT_ref));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Reference signal')

figure(3)
plot(Freq, abs(FFT_meas))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Measured signal')

figure(4)
% plot(lag, abs(corr(:,20)));
m_corr = mean(abs(corr),2);
plot(lag, m_corr);
hold on
xlabel('Cross correlation shift (GHz)')
ylabel('Amplitude (A.U.')
xlim([-20 20])

% FWHM computation (dx of points at 50% of the peak value)
[ZxP,ZxN] = find_zc(lag, m_corr,max(m_corr)/2);
FWHM = ZxN - ZxP;
plot([ZxP  ZxN], ones(2,1)*max(m_corr)/2,'r');
hold off

% some extra plots 
for ci = 1:num_windows
    m_corr = abs(corr(:,ci));
    
    figure
    hold on
    plot(lag, m_corr);
    xlabel('Cross correlation shift (GHz)')
    ylabel('Amplitude (A.U.')

    % FWHM computation (dx of points at 50% of the peak value)
    [ZxP,ZxN] = find_zc(lag, m_corr,max(m_corr)/2);
    % take only the values closest to the peak
    [pval,indp] = max(m_corr); % peak amplitude and indice
    xp = lag(indp); % x position of peak
    ind1 = find(ZxP<xp);
    ZxP = ZxP(ind1(end)); % last value
    ind2 = find(ZxN>xp);
    ZxN = ZxN(ind2(1)); % first value    
    
    FWHM = ZxN - ZxP;
    plot([ZxP  ZxN], ones(2,1)*max(m_corr)/2,'r');
    hold off

    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% 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 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
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 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end