Phase shift of two sine curves

Hello everyone,

im trying to get the phase shift of two sine curves, one is the perfect sine wave (y=sine(x)) and the other one describes a robot movement. The robot movement is sometimes faster, but often slower than the theoretical sine wave. I'm trying to figure out how to get the phase shift, but I often get the error message "Index exceeds the number of array elements (10)" because the robot movement gives me sometimes more or sometimes less then 10 y-crossings. I'm comparing the nearest numbers to the y Crossing and I think the problem is, that Matlab thinks more numbers are the y Crossing although there only can be 10. Is there a function I can use?

I tried to split the data vector into 10 parts and search for the y-Crossing with mink funktion, but then also I get two points, that are really close to each other. Here a part of my code, the robot movement data is attached. It can be assumed that the first robot data value is the first y-Crossing, it doesn't matter which number the first data value has.

z_irl_temp = z_irl;

anzahl_z_irl = numel(z_irl_temp);
teilbar = round(anzahl_z_irl * 0.1);                            %Number of elements in vector * 0.1 so I get a Matrix with 10 columns
amnachsten = anzahl_z_irl + (teilbar - rem(anzahl_z_irl,teilbar)); %This is so I can get the number of cevtor rows that are divisible 

for i = numel(z_irl):amnachsten
   
    z_irl_temp(i) = 0;                       %filling the missung elements with zeros
    
end

z_irl_temp_numel = numel(z_irl_temp);
z_irl_temp_numel = z_irl_temp_numel/10;

neu = reshape(z_irl_temp,teilbar,[]);        %reshape the new matrix to a matrix with 10 columns




const_irl = zeros(1,length(teilbar));

for i = 1:length(teilbar)
    const_irl(i) = 389.387;                     %This is where the sine wave crosses the x-axis
end

const = const';
for i = 1:width(neu)
[~,idx_irl(:,i)] = mink(abs(neu(:,i)-const_irl),1);
end

for i = 1:width(neu)
    
    if i > 1
        
        idx_irl_temp(i) = idx_irl(i) + ((i-2) * teilbar);
        
        time_yCrossings(i) = idx_irl_temp(i) * ti;
        
        %idx_
        
    else 
    
        time_yCrossings(i) = idx_irl(i) * ti;
    
    end
    
end

Here's a picture of what exactly I'm trying to find. The green line describes the ideal sine curve, the blue one the real robot data.

I'm using Matlab R2020b for academic use. Thanks for your support, I'm thankful for your help!

[dpb]

This closely mimics a Q? yesterday...see my Answer to

https://www.mathworks.com/matlabcentral/answers/868048-i-need-to-measure-the-delay-between-each-peak-of-two-signals?

Answer 1

[Mathieu NOE]

hello Bruce

the code you posted is lacking probably the section that load the data ; I could not really reproduce your plot.

so I prefered to show you my tools on dummy sinewaves - you can easily adapt it to your data.

this code generates two slighttly different sinewaves and compute the zero crossing time values for both positive and negative slope cases. Then it's fairly easy to compute the time delta between the two sets of points are shown at the end :

clc
clearvars

% dummy data
n=100;
x=(0:n-1)/n;
y1 = sin(30*x -0.5)+ 0.01*randn(1,n);
y2 = sin(31*x -0.7)+ 0.01*randn(1,n);

threshold = 0; % your value here

[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y1,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
[t0_pos2,s0_pos2,t0_neg2,s0_neg2]= crossing_V7(y2,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(1)
plot(x,y1,'b',t0_pos1,s0_pos1,'db',t0_neg1,s0_neg1,'+b','linewidth',2,'markersize',12);grid on
hold on
plot(x,y2,'g',t0_pos2,s0_pos2,'dg',t0_neg2,s0_neg2,'+g','linewidth',2,'markersize',12);grid on
hold off
legend('signal 1','signal 1 positive slope crossing points','signal 1 negative slope crossing points',...
  'signal 2','signal 2 positive slope crossing points','signal 2 negative slope crossing points'  );

% time difference for positive slope crossing points
dt_pos_slope = t0_pos1 - t0_pos2;

% time difference for negative slope crossing points
dt_neg_slope = t0_neg1 - t0_neg2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% M Noe
% added positive or negative slope condition

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            slope_sign(ii) = sign(slope);
            t0(ii) = t0(ii) - S(ind(ii)) * slope;
            s0(ii) = level;
        end
    end
end

% extract the positive slope crossing points 
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);

% extract the negative slope crossing points 
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);

end

[Bruce Rogers]

Hello Mathieu,

thanks for your answer, sorry, I just saw that I didn't say how to deal with the Data I gave you. Your solution works totally fine, but I still got the problem with the different amount of crossings with the ideal sine wave and the robot sine wave.

Matrix dimensions must agree.

Error in Frequenzdiff (line 362)
dt_neg_slope = t0_neg1 - t0_neg2;

and that's because:

How can I tell Matlab to use the first point at t=0.004 as a crossing, if the first point is not recognized? If I would always make it the first point it could happen that t0_neg1 will be a 1x6 double and I'll get the Matrix dimensions error again

[Mathieu NOE]

hello Bruce

of course my demo was simple because the number of crossing points are the same for both data sets

in your case that might be different , so you have to take only the "valid" crossing points in each data sets , for example in specifying in which valid time range they must belong

if you share your code with data I can further help you if you want

[Bruce Rogers]

Hello Mathieu,

I just used some easy if loop to check if the first crossing was detected by Matlab or not, so the whole code is running now \o/

Thanks a lot!

[Mathieu NOE]

Glad it works now !