Bi-linear (piecewise) curve fitting

26 vues (au cours des 30 derniers jours)
Mark Cuanan
Mark Cuanan le 6 Mar 2023
Can someone explain how to plot the two linear lines (Billinear curve fiiting) from the existing curve?
  20 commentaires
Mark Cuanan
Mark Cuanan le 15 Mar 2023
Déplacé(e) : John D'Errico le 30 Jan 2024
Ill get back to you on this. Thanks much
Mark Cuanan
Mark Cuanan le 2 Avr 2023
Déplacé(e) : John D'Errico le 30 Jan 2024
Updated data.. Pls help.

Connectez-vous pour commenter.

Réponse acceptée

Mathieu NOE
Mathieu NOE le 7 Mar 2023
hello
try this (a fairly simple code)
T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;
% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));
plot(x,y,'*-',x,ym,'g','Linewidth',3);
hold on
% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);
for ck = 2:idy
x2 = x(ck);
y2 = y(ck);
slope = (y2-y1)/(x2-x1);
yl = y1 + slope*(xx - x1);
% compute yellow area A
A(ck) = trapz(xx,(yy-yl));
end
% find zero crossing value for A (with linear interpolation)
threshold = 0;
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);
plot(xx,yl,'r','Linewidth',3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
  6 commentaires
Mark Cuanan
Mark Cuanan le 8 Mar 2023
as will as the coordinates on the intersection of the red and green line?
Mathieu NOE
Mathieu NOE le 13 Mar 2023
hello my friend
here you are ; code has been upgraded / expanded to answer your last questions
also I noticed I made a sign error when doing the "net" area computation. The little blue triangle I wanted to remove was with the wrong sign so instead of removing that area , I was adding it.
Now you seethe red line has a slightly steeper slope as in my previous iteration and visually you would see also that both areas are better balanced
this is the new plot as requested with intersection point coordinates and areas computation and display :
code :
T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;
% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));
xm = x(idy);
figure(1)
% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);
for ck = 2:idy
x2 = x(ck);
y2 = y(ck);
slope = (y2-y1)/(x2-x1);
yl = y1 + slope*(xx - x1);
% compute yellow area A
A(ck) = trapz(xx,(yy-yl));
% remove small triangle area above green line
% first find intersection point (lines red / green) coordinates
yint = m; % obvious
xint = (yint-y1)/slope + x1;
% plot(xint,yint,'dk',xx,yl,'r','Linewidth',1.5); % debug plot
b = xm - xint; % triangle base
h = yl(end) - m; % triangle height
At(ck) = - 0.5*b*h; % triangle area (NB : must be a negative number to be consistent with my
% definition of area A = integral of (yy-yl)
% as written above : A(ck) = trapz(xx,(yy-yl));
end
% plot A , At vs iteration index (for loop above)
figure(2)
plot(A);
hold on
plot(At);
A = A - At; % actually remove that extra triangle area (At)
plot(A);
legend('A before correction','At','A after correction');
% find zero crossing value for A (with linear interpolation)
threshold = 0;
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);
% coordinates on the intersection of the red and green line
threshold = m;
xx_intersect = find_zc(xx,yl,threshold)
yy_intersect = interp1(xx,yl,xx_intersect)
figure(1)
plot(x,y,'*-',x,ym,'g','Linewidth',2);
hold on
plot(xint,yint,'dm',xx_intersect,yy_intersect,'dk',xx,yl,'r','Linewidth',2);
text(xx_intersect,1.1*yy_intersect,['X intersect = ' num2str(xx_intersect)]);
text(xx_intersect,1.05*yy_intersect,['Y intersect = ' num2str(yy_intersect)]);
% display the area values
% compute yellow area A
Ac = cumtrapz(xx,(yy-yl));
% remove small triangle area above green line
b = xm - xx_intersect; % triangle base
h = yl(end) - m; % triangle height
At = - 0.5*b*h; % triangle area (NB : must be a negative number to be consistent with my
% definition of area A = integral of (yy-yl)
Ac = Ac - At; % actually remove that extra triangle area (At)
area_value = interp1(xx,Ac,xx_zc); % area_value is the value of Ac for x = xx_zc (point where red line crosses the data line)
text(xx_zc*0.4,y1 + slope*(xx_zc*0.5 - x1),['Area = ' num2str(area_value,'%.2e')]);
text(xx_intersect*0.97,yy_intersect*0.95,['Area = ' num2str(-area_value,'%.2e')]);
% figure(3) % debug plot for me
% plot(xx_zc,area_value,'dk',xx,Ac,'r','Linewidth',2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

Connectez-vous pour commenter.

Plus de réponses (1)

Mathieu NOE
Mathieu NOE le 3 Avr 2023
hello
here you are
all the best
for this data file the max point of the red curve is displayed on the graph
xred_max = 662.0233
yred_max = 15500;
T1 = readtable('New sample Pls open.xlsx');
x = T1.Var3;
y = T1.Var4;
% remove NaN's
idx = isnan(x);
idy = isnan(y);
id = idx & idy;
x(id) = [];
y(id) = [];
% find max force point (apex)
[m,idy] = max(y);
% ym = m*ones(size(y));
xm = x(idy);
% create the red curve for y going up to 15500
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% find 60% of max value crossing value for A (with linear interpolation)
threshold = 0.6*m;
x_c = find_zc(x,y,threshold);
y_c = interp1(x,y,x_c);
slope = (y_c-y1)/(x_c-x1);
yred_max = 15500;
xred_max = x1 + (yred_max-y1)/slope;
xred = linspace(x1,xred_max,10);
yred = y1 + slope*(xred - x1);
figure(1)
plot(x,y,'*-',xred,yred,'r','Linewidth',2);
text(xred_max*1.01,yred_max*1.01,['X = ' num2str(xred_max) ' / Y = ' num2str(yred_max)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
  2 commentaires
Mark Cuanan
Mark Cuanan le 4 Avr 2023
can u also plot the green line? Same as the previous one? Where both areas (A1 and A2) has equal values? Thanks much
Mathieu NOE
Mathieu NOE le 4 Avr 2023
hello
I wish I could be again, same remark as before , how would you solve a two unknowns problem with only one input ? the area must be equal but that is not enough to give the green line a slope and an origin
I could simply decide that the slope is zero or what is your suggestion ? how did you do in the excel file ? for me it seems you simply draw a line but how did you decide for the slope ?
all the best

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by