hello
I would suggest to compute the zero crossing time coordinates of both signals
the function used here implement linear interpolation for higher accuracy (better than one sample dt)
here a small demo of 2 signals of different amplitudes and shifted by 0.2 rad
you see the computation gives the same result as what we used in first place for the data generation
x = 0:0.1:30;
y1 = sin(x+0.25).*exp(-x/10); % reference signal
y2 = 0.75*sin(x+0.45).*exp(-x/10); % 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
