Shifting one curve to a reference curve using error minimization (fmincon from the Minimization Toolbox)

I need to shift the red curve ONLY horizontally to overlap with the black one. This should be possible using error minimization (e.g. fmincon in the Optimization Toolbox) but I have not managed to do this so far. How can I do this?

[Mathieu NOE]

hello

can you share the data ? actual data and target data ?

[Saeid]

The Data are in the excel attachment. [x1 y1] are the reference (black) and [x2 y2] are the data to be shifted horizontally. I want the shifting to be automatic, because in real life there is a large number of such datasets.

Answer 1

[Mathieu NOE]

hello again

as we look at data plotted in log log scale, a x shift is in fact a multiplicative coefficient on the Y amplitude

see the code and results below

I removed in my computation the first X samples where the two curve are not really parallel so it matches better for x > 10

clc
clearvars

data = readtable('ShiftData.xlsx');
X1 = data.x1;
Y1 = data.y1;
X2 = data.x2;
Y2 = data.y2;

% avoid taking into account first samples where two curves are not parallel
a = (Y1./Y2);
b = abs(a./rms(a));
ind = find( b> 0.5 & b < 2 );
am = rms(a(ind));

loglog(X1,Y1,X2,Y2,X2,Y2.*am);
legend('curve 1','curve 2' ,'curve 2 shifted');

[Saeid]

Hi,

thanks for the help, but rms is from the singal processing app, which I don't have. I only have Optimization and Curve Fitting. In addition, what you have done is a vertical shift, whereas what I need is a horizontal one!

Saeid

[Mathieu NOE]

oops - you're right , I was a bit too quick

now this should be better

I look for some points at same Y values belonging to both data and compute the X ratio to shift the Y2 data

here is it :

code :

clc
clearvars

data = readtable('ShiftData.xlsx');
X1 = data.x1;
Y1 = data.y1;
X2 = data.x2;
Y2 = data.y2;


ll = max(min(Y1),min(Y2));
mm = min(max(Y1),max(Y2));
ind1 = find(Y1>ll & Y1<mm);
X11 = X1(ind1);
Y11 = Y1(ind1);
ind2 = find(Y2>ll & Y2<mm);
X22 = X2(ind2);
Y22 = Y2(ind2);
ll = max(min(Y11),min(Y22));
mm = min(max(Y11),max(Y22));


N = 25;
thr = logspace(log10(3*ll),log10(0.7*mm),N);
for ci = 1: N
    [t0_pos1(ci),~,~,~]= crossing_V7(Y11,X11,thr(ci),'linear'); % positive (pos) and negative (neg) slope crossing points 
    [t0_pos2(ci),~,~,~]= crossing_V7(Y22,X22,thr(ci),'linear'); % positive (pos) and negative (neg) slope crossing points 
    dt(ci) = t0_pos2(ci)/t0_pos1(ci);

end
yy1 = interp1(X1,Y1,t0_pos1);
yy2 = interp1(X2,Y2,t0_pos2);

shift_X = 1./mean(dt);

loglog(X1,Y1,t0_pos1,yy1,'dr',X2,Y2,t0_pos2,yy2,'dc',X2.*shift_X,Y2);
legend('curve 1','curve 1 test points','curve 2' ,'curve 2 test points','curve 2 shifted');


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


% 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 ~isempty(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);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Mathieu NOE]

hello

problem solved ?

Answer 2

[Saeid]

Hi Mathieu,

thanks for the elaborate response and sorry for my late reply. Yes, it worked pretty well!

Saeid

[Mathieu NOE]

hello Saeid

would you mind accepting my answer ?

tx !

Answer 3

[Saeid]

Sure. I just did!