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

2 次查看（过去 30 天）
Saeid
Saeid 2021-11-5
回答： Saeid 2021-11-22
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?
  2 个评论
Saeid
Saeid 2021-11-5
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Mathieu NOE
Mathieu NOE 2021-11-5
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');
  3 个评论
显示 1更早的评论
Mathieu NOE
Mathieu NOE 2021-11-5
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

请先登录，再进行评论。

更多回答（2 个）

Saeid
Saeid 2021-11-21
Hi Mathieu,
thanks for the elaborate response and sorry for my late reply. Yes, it worked pretty well!
Saeid
Saeid
Saeid 2021-11-22
Sure. I just did!

类别

在 Help Center 和 File Exchange 中查找有关 Least Squares 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by